cgcm.cpp

数学计算程序

		}
		else if(strcmp(argv[i], "-outpdb") == 0){
			i++;
			if(i >= argc){
				cerr << argv[0] << ": -outpdb needs a filename.\n";
				exit(1);
			}
			else{
				str << argv[i];
				str >> outpdb;
			}
		}
		else if(strcmp(argv[i], "-outtxt") == 0){
			i++;
			if(i >= argc){
				cerr << argv[0] << ": -outtxt needs a filename.\n";
				exit(1);
			}
			else{
				str << argv[i];
				str >> outtxt;
			}
		}
		else if(strcmp(argv[i], "-center") == 0){
			i++;
			if(i >= argc){
				cerr << argv[0] << ": -outtxt needs a filename.\n";
				exit(1);
			}
			else{
				str << argv[i];
				str >> center;
			}
		}
		else if(strcmp(argv[i], "-element") == 0){
			i++;
			if(i >= argc){
				cerr << argv[0] << ": -element needs a filename.\n";
				exit(1);
			}
			else{
				str << argv[i];
				str >> element;
			}
		}
		else if(strcmp(argv[i], "-parafile") == 0){
			i++;
		}
		else{
			cerr << argv[0] << ": unrecognized arguement:"
				<< argv[i] << "\n\n";
			cerr << "parameters include: \n"
				 << " -grad_tol  -g_rho_min  -g_rho_max  -g_phi_min  -g_phi_max \n"
				 << " -g_z_min   -g_z_max    -rho_scale  -z_scale    -center \n"
				 << " -outpdb    -outtxt     -element    -parafile \n";
			exit(1);
		}
	}

	////////////////////////////////////////////////////////////
	/// read .pdb file to obtain number of atoms
	num = pdb_num(infile);
	dim = num*3;
	atoms = new CParticles(num);
	////////////////////////////////////////////////////////////

	ofstream fout(outtxt.c_str());
	if(!fout)
	{
		cerr << "-------------------------------------------------------" << endl;
		cerr << "Warning: Cannot write to " << outtxt.c_str()
			 << " for output." << endl;
		cerr << "-------------------------------------------------------" << endl;
	}

	double *x = new double[num];
	double *y = new double[num];
	double *z = new double[num];
	double *low = new double[dim];
	double *up = new double[dim];
	double *var = new double[dim];

	/// change epsilon  Parm->eps = 1.e-8;
	/// line 698 file asa_cg.cpp
	/*asacg_parm cgParm ;
	asa_parm asaParm ;
	asa_cg_default (&cgParm) ;
	asa_default (&asaParm) ;
	cgParm.eps = epsilon;
	asaParm.eps = epsilon;*/

	read_pdb(infile, x, y, z, element);
	if(center == 1) center_cyl(x, y, z, num);
	//write_pdb(string("testtest.pdb"), x, y, z, num, string("Ni"));

	atoms->Read(x, y, z);
	atoms->WriteCyl(x, y, z);

	_3d_1d(x, y, z, num, var);
	/*_1d_3d(var, x, y, z, num);
	write_pdb(string("cyl_xyz2.pdb"), x, y, z, num, string("Ni"));*/

	/// Boundary constraints
	atoms->SetBound(low, up, g_rho_min, g_rho_max, g_phi_min, g_phi_max, g_z_min, g_z_max);
	/// parameter corresponding to element
	setup_param(element);
	for(int i = 0; i < dim;)
	{
		var[i] *= rho_scale;//0.95;
		var[i+2] *= z_scale;
		i += 3;
	}

	///////////////////////////////////////////////////////////////////////////
	/// simple boundary test module, you need to improve it later.
	double rhom, zm;
	boundary_cyl(var, dim, rhom, zm);
	if(rhom > g_rho_max) 
	{
		cerr << "-------------------------------------------------------" << endl;
		cerr << "Warning: rho_max = " << rhom 
			 << " of atoms out of constraints " << g_rho_max << endl;
		cerr << "Suggestion: set rho_scale smaller than " << g_rho_max/rhom << endl;
		cerr << "      But the algorithm may not be tolerant. " << endl;
		cerr << "-------------------------------------------------------" << endl;
	}
	if(zm > g_z_max) 
	{
		cerr << "-------------------------------------------------------" << endl;
		cerr << "Warning: z_max " << zm 
			 << "of atoms out of constraints " << g_z_max << endl;
		cerr << "Suggestion: set z_scale smaller than " << g_z_max/zm << endl;
		cerr << "      But the algorithm may not be tolerant. " << endl;
		cerr << "-------------------------------------------------------" << endl;
	} 
	/// simple boundary, you need to improve it later
	///////////////////////////////////////////////////////////////////////////

	//asa_cg (var, low, up, dim, NULL, NULL, NULL, grad_tol, EnergySC_CYL, ForceSC_CYL, ForceEnergySC_CYL, NULL);

	asa_cg (var, low, up, dim, NULL, NULL, NULL, grad_tol, EnergySC_CYL, ForceSC_CYL, NULL, NULL);

	_1d_3d(var, x, y, z, num);
	atoms->ReadCyl(x, y, z);
	atoms->Write(x, y, z);

	//// //// Test if energy caculated by EnergySC and ForceSC is equal.
	//int flag = 0; //////////
	//for(int i = 0; i < num; i++)/////////
	//{////////
	//	if(fabs(x[i]-x1[i]) > 1.e-3) flag = 1;///
	//	if(fabs(y[i]-y1[i]) > 1.e-3) flag = 1;///
	//	if(fabs(z[i]-z1[i]) > 1.e-3) flag = 1;///
	//}////////
	//cout << "flag: " << flag << endl;//////
	//cout << "energy: " << EnergySC_CYL(var, dim) << endl;/////////////
	//ForceSC_CYL(ff, var, dim);/////////////
	////for(int i = 0; i < dim; i ++)///////////
	//	//cout << ff[i] << "  ";/////////////
	//cout << "ff energy: " << ForceSC(x, y, z, num, fff, string("Ni")) << endl;//////////

	write_pdb(outpdb, x, y, z, num, element);
	cout << num << " atoms" << endl << endl;
	fout << num << " atoms" << endl << endl;

	fout.close();
	delete []x;
	delete []y;
	delete []z;
	delete []var;
	delete []low;
	delete []up;
	delete atoms;

	//delete []x1;/////
	//delete []y1;////
	//delete []z1;////
	//delete []ff;////
	//for(int i = 0; i < 3; i++) delete []fff[i];/////
	//delete []fff;//////

	//double xxx[2] = {0,0};//////////////////////////////////////////////////////////////////////////
	//double yxx[2] = {0,0};
	//double zxx[2] = {0, 0.09 };
	//double **ffff = new double*[3];
	//for(int i=0; i<3; i++) ffff[i] = new double[2];
	//
	//cout << EnergySC(xxx, yxx, zxx, 2, string("Ni"))<<endl;
	//ForceSC(xxx,yxx,zxx,2,ffff,string("Ni"));
	//for(int i=0; i < 3; i++)
	//{//////////////////////////////////////////////////////////////////////////
	//	for(int j=0; j<2 ; j++)
	//	{
	//		cout << ffff[i][j] << "  ";
	//	}
	//	cout << endl;
	//}
	//for(int i=0; i<3; i++) delete []ffff[i];
	//delete []ffff;//////////////////////////////////////////////////////////////////////////
}

/* evaluate the objective function */
double EnergySC_CYL(double *xx, int n)
{
    double f;
    double *x = new double[n/3];
	double *y = new double[n/3];
	double *z = new double[n/3];
	for(int i = 0; i < n;)
	{
		x[i/3] = xx[i];
		y[i/3] = xx[i+1];
		z[i/3] = xx[i+2];
		i += 3;
	}
	//setup_param(string("Ni"));
	/// cylindrical coordination to cartesian coordination
	atoms->ReadCyl(x, y, z);
	atoms->Write(x, y, z);
	//Rc = 100;
	f = EnergySC(x, y, z, n/3, string("Ni"));
	delete []x;
	delete []y;
	delete []z;
    return (f);
}

/* evaluate the gradient of the objective function */
void ForceSC_CYL(double  *g, double  *xx, int n)
{
	int i;
	double phi;
    double *x = new double[n/3];
	double *y = new double[n/3];
	double *z = new double[n/3];
	double **f = new double*[3];
	for(i = 0; i < 3; i++) f[i] = new double[n/3];
	for(i = 0; i < n;)
	{
		x[i/3] = xx[i];
		y[i/3] = xx[i+1];
		z[i/3] = xx[i+2];
		i += 3;
	}
	//setup_param(string("Ni"));
	//Rc = 100;
	/// cylindrical coordination to cartesian coordination
	/// need converted matrix
	atoms->ReadCyl(x, y, z);
	atoms->Write(x, y, z);
	ForceSC(x, y, z, n/3, f, string("Ni"));
	for(i = 0; i < n/3; i++)
	{
		phi = atoms->atoms[i]->cyl.phi;
		g[3*i] = -( f[0][i]*cos(phi) + f[1][i]*sin(phi) );
		g[3*i+1] = -( -f[0][i]*sin(phi) + f[1][i]*cos(phi) );
		g[3*i+2] = -( f[2][i] );
	}
	delete []y;
	delete []z;
	delete []x;
	for(i = 0; i < 3; i++) delete []f[i];
	delete []f;
}

/* value and gradient of the objective function */
double ForceEnergySC_CYL (double  *g, double  *xx, INT32 n)
{
	//double energy = ForceSC_CYL(double  *g, double  *xx, int n);
	int i;
	double energy, phi;
    double *x = new double[n/3];
	double *y = new double[n/3];
	double *z = new double[n/3];
	double **f = new double*[3];
	for(i = 0; i < 3; i++) f[i] = new double[n/3];
	for(i = 0; i < n;)
	{
		x[i/3] = xx[i];
		y[i/3] = xx[i+1];
		z[i/3] = xx[i+2];
		i += 3;
	}
	//setup_param(string("Ni"));
	//Rc = 100;
	/// cylindrical coordination to cartesian coordination
	atoms->ReadCyl(x, y, z);
	atoms->Write(x, y, z);
	energy = ForceSC(x, y, z, n/3, f, string("Ni"));
	for(i = 0; i < n/3; i++)
	{
		phi = atoms->atoms[i]->cyl.phi;
		g[3*i] = -( f[0][i]*cos(phi) + f[1][i]*sin(phi) );
		g[3*i+1] = -( -f[0][i]*sin(phi) + f[1][i]*cos(phi) );
		g[3*i+2] = -( f[2][i] );
	}
	delete []y;
	delete []z;
	delete []x;
	for(i = 0; i < 3; i++) delete []f[i];
	delete []f;
	return energy;
}