C++ Molecule::DistanceOfMassCenter方法代码示例

double NWenergy(Molecule a, Molecule b, string basis = "6-31G")
{
	string filename("DATA/NW.nw");
	ToNWchemFileHF(a,b,filename,basis);
	//Use shell script to solve the scf energy
	double total_energy = 0;
#ifdef _WIN32
	cout << "It is tested under win32 os." << endl;
	double r = a.DistanceOfMassCenter(b);
	total_energy = 33*(1/r/r/r/r/r/-1/r/r);
#else
	system("./NW_perform.sh");
	total_energy = ReadFile("DATA/temp.txt");
#endif
	CallTimes += 1;
	return total_energy;
}

void MonteCarlo(Molecule a, Molecule b, double StepLength, double StepPrecision, double MaxRotTime, double RotPrecision, string forcefieldORbasis, string MinConfigName,string RelativeMinConfigName)
{
	//initi parameter
	//StepPrecision = 0.15;
	//StepLength = 6;
	const int stepCountTimes = (int)(StepLength / StepPrecision) + 1;
	//RotPrecision = 10;
	//MaxRotTime = 10000;
	const double OriginTemp=400;
	const double FirstDistance = 15;
	double dT=7;
	double Temp = OriginTemp;
        int ControlNum=0;
        const int MaxControlNum=(int)((360/RotPrecision)*(360/RotPrecision)*10);

	int MinConfigNum = 0;
	int RelativeMinConfigNum = 0;

	outputWelcome();
	//Calculate zeroEnergy
	//const double ZeroEnergy = CalculateEnergy(a, forcefieldORbasis) + CalculateEnergy(b, forcefieldORbasis);
	const double ZeroEnergy = 0;


	/*initilize min value with first step opt*/
	//Initilized 2 molecules
	int firstStepOptTimes = (int)((FirstDistance-4)/StepPrecision);
	Eigen::Vector3d myFirstTrans;
	myFirstTrans << FirstDistance, 0, 0;// set 2 molecules as a far distance in the begining
	InitiConfig(a, b, myFirstTrans);
	//Calculate initi potential
	double potential = CalculatePotential2(a,b,forcefieldORbasis,ZeroEnergy);

	//Initilize relativeMIn, relativeMIN is for each step-opt procedure
	double relativeMinPotential = potential;
	DoubleMolecule relativeMinConfig(a,b,relativeMinPotential);
	//begin move in step-opt
	for (int i = 0; i != firstStepOptTimes; i++)
	{
		b.PerformXTrans(-1*StepPrecision);
		potential = CalculatePotential2(a, b, forcefieldORbasis, ZeroEnergy);
		if (potential < relativeMinPotential)
		{
			relativeMinPotential = potential;
			relativeMinConfig.Set(a,b,relativeMinPotential);
		}
	}
	//Initilize global Min config & potential
	double MinPotential = relativeMinPotential;
	DoubleMolecule MinConfig(relativeMinConfig);
	MinConfig.ToXYZ(CombineFileName(RelativeMinConfigName, RelativeMinConfigNum), CallTimes);
	RelativeMinConfigNum += 1; MinConfigNum += 1;//Save this one
	cout << "No." << MinConfigNum << " is obtained " << endl;
	MinConfig.output();



	/*Random rot*/

	for (int i = 0; i != MaxRotTime; i++)
	{
                //Judge: 
               if(a.DistanceOfMassCenter(b)>FirstDistance)
                      break;
	        Temp+=dT; ControlNum+=1;
                if(ControlNum>MaxControlNum)
                           break;
		//Set a,b as step-opt. This is important for each loop
		MinConfig.GetInfo(a, b, potential);//Form global min config, to perform rot
		PerformRandomRotEuler2(a, b, RotPrecision);
		potential = CalculatePotential2(a, b, forcefieldORbasis, ZeroEnergy);
		if (MonteCarlo01Distribution(potential - MinPotential, Temp) !=0)
		{
		    cout<<"Random rot from global min configuration is permitted!, T is "<<Temp<<endl;
			//step-opt
			StepOpt(a, b, relativeMinPotential, relativeMinConfig, potential, ZeroEnergy, forcefieldORbasis, StepLength, StepPrecision);
			relativeMinConfig.ToXYZ(CombineFileName(RelativeMinConfigName, RelativeMinConfigNum),CallTimes); RelativeMinConfigNum += 1;//Save this relative one
			cout << "No." << RelativeMinConfigNum << " Relative MinConfig rot from No." << MinConfigNum<<" global min config and potential is " << relativeMinPotential << endl;
			//If after-rot step-opt config has lower energy, save as global min
			if (relativeMinPotential < MinPotential)
			{
				MinPotential = relativeMinPotential;
				MinConfig = relativeMinConfig;
				Temp=OriginTemp;//After obtain a global min config, set T as origin temperature
				cout << "No." << MinConfigNum << " MinConfig obtained from Relative MinConfig" << endl;
				MinConfig.ToXYZ(CombineFileName(MinConfigName, MinConfigNum), CallTimes); MinConfigNum += 1;//Save this one
				MinConfig.output();
                                 ControlNum=0;
			}
			//If the relative min is not global min, keep rot from relative min until forbidden!
			else
			{
                double BranchTemp=Temp;
				for (int iBranch = 0; iBranch != MaxRotTime; iBranch++)
				{
					relativeMinConfig.GetInfo(a, b, potential); relativeMinPotential = potential;//set a,b as step-opt config from relative min
					//perform random rot
					//BranchTemp initilize as Temp, and begin to cool down when rot from relative min
                                        if(BranchTemp>0.8*OriginTemp) BranchTemp-=2*dT;
					PerformRandomRotEuler2(a, b, RotPrecision);
//.........这里部分代码省略.........