C++ CMol::Cal_E方法代码示例

void Cal_E_MM_Scaned(void)
{
	int i, j, nAtom, Count=0;
	double E_Min_Local_QM;

	E_Scan_MM_Mean = 0.0;
	E_Scan_QM_Mean = 0.0;
	nAtom = Mol_ESP.nAtom;

	for(i=0; i<n_Phi; i++)	{
		E_Min_Local_QM = 1.0E100;
		for(j=0; j<nScan_List[i]; j++)	{
			if(E_Scan_QM[i][j] < E_Min_Local_QM)	{	// to find the lowest QM energy
				E_Min_Local_QM = E_Scan_QM[i][j];
			}
		}

		for(j=0; j<nScan_List[i]; j++)	{
			memcpy(Mol_ESP.x, x_Scan_list[i][j], sizeof(double)*nAtom);
			memcpy(Mol_ESP.y, y_Scan_list[i][j], sizeof(double)*nAtom);
			memcpy(Mol_ESP.z, z_Scan_list[i][j], sizeof(double)*nAtom);
			E_Scan_MM[i][j] = Mol_ESP.Cal_E(0);

			if( E_Scan_QM[i][j] < (E_Min_Local_QM+6.0) )	{	// skip configuration with much higher energies
				E_Scan_QM_Mean += E_Scan_QM[i][j];
				E_Scan_MM_Mean += E_Scan_MM[i][j];
				Count ++;
			}
		}
	}
	E_Scan_QM_Mean /= Count;
	E_Scan_MM_Mean /= Count;
}

double func_Geo_Opt_Constraint(unsigned n, const double *x, double *grad, void *my_func_data)
{
	int i, iPos, nAtom;
	double E_Total, *gx, *gy, *gz, *xx, *yy, *zz;
	CMol* pMol;

	pMol = (CMol*)my_func_data;
	nAtom = pMol->nAtom;
	gx = pMol->grad_x;
	gy = pMol->grad_y;
	gz = pMol->grad_z;
	xx = pMol->x;
	yy = pMol->y;
	zz = pMol->z;

	if(grad)	{
		Iter_Opt_Constrained ++;
		for(i=0; i<nAtom; i++)	{	// update the coordinate
			iPos = 3*i;
			xx[i] = x[iPos  ];
			yy[i] = x[iPos+1];
			zz[i] = x[iPos+2];
		}

		E_Total = pMol->Cal_E(0);

		for(i=0; i<nAtom; i++)	{	// get the gradient
			iPos = 3*i;
			grad[iPos  ] = gx[i];
			grad[iPos+1] = gy[i];
			grad[iPos+2] = gz[i];
		}
	}
	else	{
		E_Total = pMol->Cal_E(0);
	}

//  printf( "TT: %f\n", E_Total );
	return E_Total;
}

void Output_1D_QM_MM(void)
{
	FILE *fOut;
	char szName[256];
	int i, j, nAtom;
	double E_Min_1D_QM, E_Min_1D_MM_Mean, E_Min_1D_QM_Mean, E_Min_1D_QM_MM_Shift, Weight_Sum;

	if(ProgID != 0)	{
		return;
	}

	nAtom = Mol_ESP.nAtom;
	for(i=0; i<n_Phi; i++)	{
		sprintf(szName, "1d-qm-mm-%d.dat", i+1);
		fOut = fopen(szName, "w");

		E_Min_1D_QM = 1.0E100;
		for(j=0; j<nScan_List[i]; j++)	{
			memcpy(Mol_ESP.x, x_Scan_list[i][j], sizeof(double)*nAtom);
			memcpy(Mol_ESP.y, y_Scan_list[i][j], sizeof(double)*nAtom);
			memcpy(Mol_ESP.z, z_Scan_list[i][j], sizeof(double)*nAtom);
			E_Scan_MM[i][j] = Mol_ESP.Cal_E(0);

			if(E_Min_1D_QM > E_Scan_QM[i][j])	{
				E_Min_1D_QM = E_Scan_QM[i][j];
			}
		}

		E_Min_1D_MM_Mean = E_Min_1D_QM_Mean = 0.0;
		Weight_Sum = 0.0;
		for(j=0; j<nScan_List[i]; j++)	{
			Weight_Sum += w_Scan_List[i][j];
			E_Min_1D_MM_Mean += (w_Scan_List[i][j] * E_Scan_MM[i][j] );
			E_Min_1D_QM_Mean += (w_Scan_List[i][j] * E_Scan_QM[i][j] );
		}
		E_Min_1D_QM_MM_Shift = (E_Min_1D_QM_Mean - E_Min_1D_MM_Mean)/Weight_Sum;


		for(j=0; j<nScan_List[i]; j++)	{
			fprintf(fOut, "%.1lf %.5lf %.5lf\n", Phi_Scan_List[i][j], E_Scan_QM[i][j]-E_Min_1D_QM, E_Scan_MM[i][j]+E_Min_1D_QM_MM_Shift-E_Min_1D_QM);	// the relative 1D energy
		}

		fclose(fOut);
	}
}

void MM_Fixed_1D_Scan(void)
{
	int i, Phi, Count;

	RestoreCoordinates();
	for(i=0; i<n_Phi; i++)	{
		Mol.QueryDihedral(IdxDihSelect[i]);
		Mol.Edit_Dihedral(IdxDihSelect[i], Phi_To_Set[i]);
	}
	
	Count = 0;
	for(Phi=-180; Phi<=180; Phi+=BIN_SIZE)	{	// 1D scan for this rotamer
		Phi_To_Set[Active_Phi] = Phi*1.0;
		Mol.QueryDihedral(IdxDihSelect[Active_Phi]);
		Mol.Edit_Dihedral(IdxDihSelect[Active_Phi], Phi_To_Set[Active_Phi]);
		E_Phi[Count] = Mol.Cal_E(0);
		Count++;
	}

	return;
}

void Extract_Coord_E(char szName[], int Idx_Phi, int Idx_Part)
{
	FILE *fIn, *fOut;
	int i, nAtom, ToRead=1, ReadItem, n_Rec=0;
	char szOutput[256], szErrorMsg[256], szLine[256], *ReadLine, szTag[256], ErrorMsg[256];
	double E_QM, E_MM, Phi, Phi_Set, x_Save[MAX_ATOM], y_Save[MAX_ATOM], z_Save[MAX_ATOM];

	nAtom = Mol.nAtom;
	memcpy(x_Save, Mol.x, sizeof(double)*nAtom);
	memcpy(y_Save, Mol.y, sizeof(double)*nAtom);
	memcpy(z_Save, Mol.z, sizeof(double)*nAtom);


	sprintf(szOutput, "tor-1D-idx-%d.dat", Idx_Phi+1);
	fOut = fopen(szOutput, "a+");
	fseek(fOut, 0, SEEK_END);

	fIn = fopen(szName, "r");
	if(fIn == NULL)	{
		sprintf(szErrorMsg, "Fail to open %s\nQuit\n", szName);
		Quit_With_Error_Msg(szErrorMsg);
	}


	while(ToRead)	{
		if(feof(fIn))	{
			break;
		}
		ReadLine = fgets(szLine, 256, fIn);
		if(ReadLine == NULL)	{
			break;
		}
		else	{
//			if(FindString(szLine, " Center     Atomic      Atomic")>=0)	{	// to extract the coordinate
//				Skip_N_Line(fIn, szName, 2);
//			if(FindString(szLine, "     Input orientation:")>=0)	{	// to extract the coordinate
			if(FindString(szLine, " orientation:")>=0)	{	// to extract the coordinate
				Skip_N_Line(fIn, szName, 4);
				for(i=0; i<nAtom; i++)	{
					ReadLine = fgets(szLine, 256, fIn);
					if(ReadLine == NULL)	{
						break;
					}
					ReadItem = sscanf(szLine+31, "%lf %lf %lf", &(Mol.x[i]), &(Mol.y[i]), &(Mol.z[i]));
					if(ReadItem != 3)	{
						ToRead = 0;
						break;
					}
				}
			}
			else if( (FindString(szLine, " SCF Done: ")>=0) && (QM_Level == QM_LEVEL_HF) )	{	// HF
				E_QM = Get_Energy(szLine);
			}
			else if( (FindString(szLine, "EUMP2 =")>=0) && (QM_Level == QM_LEVEL_MP2) )	{
				E_QM = Get_Energy(szLine+27);
			}
			else if(FindString(szLine, " Optimization completed")>=0)	{
				sprintf(szTag, "  D(%d,%d,%d,%d)", DihList[Idx_Phi][0]+1, DihList[Idx_Phi][1]+1, DihList[Idx_Phi][2]+1, DihList[Idx_Phi][3]+1);
				To_Find_Tag(fIn, szName, szTag, szLine);
				ReadItem = sscanf(szLine+28, "%lf", &Phi_Set);	// previous, modredundant
				if(ReadItem != 1)	{
					sprintf(ErrorMsg, "Error in extracting the dihedral.\n%s\nQuit\n", szLine);
					Quit_With_Error_Msg(ErrorMsg);
				}

				Phi = Mol.QueryDihedral(IdxDihSelect[Idx_Phi]);
				E_MM = Mol.Cal_E(0);

				fprintf(fOut, "E_Scan %.13E %.13E  Phi  %.1lf\n", E_QM, E_MM, Phi);
				fprintf(fOut, "Coordinate\n");
				for(i=0; i<nAtom; i++)	{
					fprintf(fOut, "%12.6lf %12.6lf %12.6lf\n", Mol.x[i], Mol.y[i], Mol.z[i]);
				}

				n_Rec++;
			}
		}
	}




	fclose(fOut);
	fclose(fIn);

	memcpy(Mol.x, x_Save, sizeof(double)*nAtom);
	memcpy(Mol.y, y_Save, sizeof(double)*nAtom);
	memcpy(Mol.z, z_Save, sizeof(double)*nAtom);
}

double CalObjectiveFunction(double Para_List[], int Get_E_Rotamer)
{
	int i, j, nAtom, time_1, time_2;
	double d_E, E_Shift, w_1D;
	double Chi_SQ_Rotamer_Local;
	double E_Rotamer_MM_Local[MAX_ROTAMER];
	double E_Min_1D_MM_Mean, E_Min_1D_QM_Mean, E_Min_1D_QM_MM_Shift, Weight_Sum;

	Chi_SQ = Chi_SQ_1D = Chi_SQ_Rotamer = 0.0;
	Chi_SQ_Rotamer_Local = 0.0;

	w_1D = 1.0 - w_Rotamer;

	time_1 = time(NULL);

	Distribute_Torsion_Parameters();
	Assign_Torsion_Parameters();

	E_Shift = Para_List[10*n_Phi];

	nAtom = Mol_ESP.nAtom;
	for(i=0; i<n_Phi; i++)	{
		for(j=0; j<nScan_List[i]; j++)	{
			memcpy(Mol_ESP.x, x_Scan_list[i][j], sizeof(double)*nAtom);
			memcpy(Mol_ESP.y, y_Scan_list[i][j], sizeof(double)*nAtom);
			memcpy(Mol_ESP.z, z_Scan_list[i][j], sizeof(double)*nAtom);
			E_Scan_MM[i][j] = Mol_ESP.Cal_E(0);
		}

		E_Min_1D_MM_Mean = E_Min_1D_QM_Mean = 0.0;
		Weight_Sum = 0.0;
		for(j=0; j<nScan_List[i]; j++)	{
			Weight_Sum += w_Scan_List[i][j];
			E_Min_1D_MM_Mean += (w_Scan_List[i][j] * E_Scan_MM[i][j] );
			E_Min_1D_QM_Mean += (w_Scan_List[i][j] * E_Scan_QM[i][j] );
		}
		E_Min_1D_QM_MM_Shift = (E_Min_1D_QM_Mean - E_Min_1D_MM_Mean)/Weight_Sum;

		for(j=0; j<nScan_List[i]; j++)	{
			d_E = (E_Scan_QM[i][j] - E_Scan_MM[i][j] - E_Min_1D_QM_MM_Shift);	// fit the relative 1D energy profile
			Chi_SQ_1D += (d_E*d_E*w_Scan_List[i][j]);
		}
	}

	
	Chi_SQ_1D *= (w_1D/w_Scan_Sum);


	for(i=0; i<nRotamer; i++)	{
		E_Rotamer_MM[i] = E_Rotamer_MM_Local[i] = 0.0;
		if(i%nProc != ProgID)	{	//only proceed when i%nProc == ProgID, do job decomposition
			continue;
		}
		memcpy(Mol_ESP.x, x_Rotamer_Opt[i], sizeof(double)*nAtom);
		memcpy(Mol_ESP.y, y_Rotamer_Opt[i], sizeof(double)*nAtom);
		memcpy(Mol_ESP.z, z_Rotamer_Opt[i], sizeof(double)*nAtom);
		Mol_ESP.FullGeometryOptimization_LBFGS_step(1);
		E_Rotamer_MM_Local[i] = Mol_ESP.Cal_E(0);
		d_E = E_Rotamer_QM[i] - E_Rotamer_MM_Local[i] - E_Shift;
		Chi_SQ_Rotamer_Local += (d_E*d_E*w_rotamer_List[i]);
	}
	if(Get_E_Rotamer)	{
//eliot (from Matt Harvey)
	for( int i = 0 ; i < nRotamer; i++ ) { E_Rotamer_MM[i] += E_Rotamer_MM_Local[i]; }
                Chi_SQ_Rotamer = Chi_SQ_Rotamer_Local;

//eliot		MPI_Allreduce(E_Rotamer_MM_Local, E_Rotamer_MM, nRotamer, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
	}
//eliot	MPI_Allreduce(&Chi_SQ_Rotamer_Local, &Chi_SQ_Rotamer, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
//	MPI_Bcast(&Chi_SQ_Rotamer, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

	printf(" w_rotamer_Sum = %f\n", w_rotamer_Sum );

	Chi_SQ_Rotamer *= (w_Rotamer/w_rotamer_Sum);

	Chi_SQ = Chi_SQ_1D + Chi_SQ_Rotamer;

//	MPI_Barrier(MPI_COMM_WORLD);
//	printf("%lf\n", Chi_SQ);
//	MPI_Finalize();

	time_2 = time(NULL);
/*
	if(ProgID == 0)	{
		FILE *fOut;
		fOut = fopen("time-rec.txt", "a+");
		fseek(fOut, 0, SEEK_END);
		fprintf(fOut, "%d\n", time_2-time_1);
		fclose(fOut);
	}
*/
	printf("Chi_SQ=%f\n", Chi_SQ );
	return Chi_SQ;
}

void Cal_E_MM_QM_Diff(int Idx)
{
	FILE *fIn, *fOut;
	char szName[256], szLine[256], *ReadLine, szTmp[256], ErrorMsg[256];
	double E_QM_Min, E_MM_Min, E_MM_Read, E_MM_Cal;
	double Para_k_Dih_Save[6], Para_phi_Dih_Save[6];
	double x_Save[MAX_NUM_ATOM], y_Save[MAX_NUM_ATOM], z_Save[MAX_NUM_ATOM];
	int i, j, ReadItem, nAtom;

	memcpy(Para_k_Dih_Save, Mol_ESP.Para_k_Dih[IdxDihSelect[Idx]]+1, sizeof(double)*6);
	memcpy(Para_phi_Dih_Save, Mol_ESP.Para_phi[IdxDihSelect[Idx]]+1, sizeof(double)*6);
	memset(Mol_ESP.Para_k_Dih[IdxDihSelect[Idx]]+1, 0, sizeof(double)*6);	// turn of this dihedral
	memset(Mol_ESP.Para_phi[IdxDihSelect[Idx]]+1, 0, sizeof(double)*6);

	sprintf(szName, "tor-1D-idx-%d.dat", Idx+1);
	fIn = fopen(szName, "r");
	if(fIn == NULL)	{
		sprintf(ErrorMsg, "Fail to open file: %s\nQuit\n", szName);
		printf("%s", ErrorMsg);
		Quit_With_Error_Msg(ErrorMsg);
	}

	sprintf(szName, "mm-tor-1D-idx-%d.dat", Idx+1);	// the one with optimized geometry in MM force field
	fOut = fopen(szName, "w");

	nAtom = Mol_ESP.nAtom;
	nScan = 0;
	while(1)	{
		if(feof(fIn))	{
			break;
		}
		ReadLine = fgets(szLine, 256, fIn);
		if(ReadLine == NULL)	{
			break;
		}
		if(strncmp(szLine, "E_Scan", 6)==0)	{
			ReadItem = sscanf(szLine, "%s %lf %lf %s %lf", szTmp, &(E_QM_Scaned[nScan]), &E_MM_Read, szTmp, &(Phi_Scaned[nScan]));


			if(ReadItem != 5)	{
				sprintf(ErrorMsg, "Error in reading file: %s\nQuit\n", szName);
				Quit_With_Error_Msg(ErrorMsg);
			}

			ReadLine = fgets(szLine, 256, fIn);	// skip one line, "Coordinate"

			for(i=0; i<nAtom; i++)	{	// read one snapshot
				ReadLine = fgets(szLine, 256, fIn);
				ReadItem = sscanf(szLine, "%lf %lf %lf", &(Mol_ESP.x[i]), &(Mol_ESP.y[i]), &(Mol_ESP.z[i]));
			}
			E_MM_Cal = Mol_ESP.Cal_E(0);	// the nergy without optimization and dihedral term is turned off

			To_Setup_All_Dihedral_Constraint();
			E_MM_Scaned[nScan] = Geoometry_Optimization_With_Constraint(&Mol_ESP);	// do restrained geometry optimizatio. soft restraint on all dihedrals. Rigid restraints on soft dihedrals

			if(E_MM_Scaned[nScan] > (1.0E3+E_MM_Cal))	{	// not a valid configuration. Something wrong in the constrained optimization
        printf(" Bad configuration: %d %f %f\n", nScan, E_MM_Scaned[nScan], E_MM_Cal );
				continue;
			}

			memcpy(x_Save, Mol_ESP.x, sizeof(double)*nAtom);
			memcpy(y_Save, Mol_ESP.y, sizeof(double)*nAtom);
			memcpy(z_Save, Mol_ESP.z, sizeof(double)*nAtom);
			
			fprintf(fOut, "E_Scan %.13E %.13E Phi %.1lf\n", E_QM_Scaned[nScan], E_MM_Scaned[nScan], Phi_Scaned[nScan]);
			fprintf(fOut, "Coordinate\n");
			for(j=0; j<nAtom; j++)	{
				fprintf(fOut, "%14.6lf%14.6lf%14.6lf\n", x_Save[j], y_Save[j], z_Save[j]);
			}
			//		fprintf(fOut, "%.1lf %.6lf %.6lf %.6lf\n", Phi_Scaned[i], E_QM_Scaned[i], E_MM_Scaned[i], E_QM_Scaned[i] - E_MM_Scaned[i]);
			

			nScan++;
		}

	}


	fclose(fOut);

	fclose(fIn);

	E_QM_Min = E_MM_Min = 1.0E100;
	for(i=0; i<nScan; i++)	{
		if(E_QM_Min > E_QM_Scaned[i])	{
			E_QM_Min = E_QM_Scaned[i];
		}
		if(E_MM_Min > E_MM_Scaned[i])	{
			E_MM_Min = E_MM_Scaned[i];
		}
	}

	for(i=0; i<nScan; i++)	{
		E_QM_Scaned[i] = (E_QM_Scaned[i]-E_QM_Min)*HARTREE_To_KCAL;
		E_MM_Scaned[i] -= E_MM_Min;
		E_QM_MM_Diff[i] = E_QM_Scaned[i]-E_MM_Scaned[i];

		if(E_MM_Scaned[i] > 60.0)	{	// an unusual large MM energy
			w_Point[i] = 0.0;
		}
//.........这里部分代码省略.........