C++ CMol类代码示例

double Geo_Constraint_Bond(unsigned n, const double *x, double *grad, void *data)
{
	int ia, ib, iPos;
	double r, r0, g[2][3];
	CMol* pMol;
	GEO_FIX_BOND *Geo_Fix_r;

	Geo_Fix_r = (GEO_FIX_BOND*)data;
	pMol = (CMol*)(Geo_Fix_r->pMol);
	ia = Geo_Fix_r->ia;
	ib = Geo_Fix_r->ib;
	r0 = Geo_Fix_r->r0;

	if(grad)	{
		r = pMol->Query_Distance(ia, ib, 1, g);
		memset(grad, 0, sizeof(double)*n);

		iPos = 3 * ia;
		grad[iPos  ] = g[0][0];
		grad[iPos+1] = g[0][1];
		grad[iPos+2] = g[0][2];

		iPos = 3 * ib;
		grad[iPos  ] = g[1][0];
		grad[iPos+1] = g[1][1];
		grad[iPos+2] = g[1][2];
	}
	else	{
		r = pMol->Query_Distance(ia, ib, 0, g);
	}

	return (r - r0);
}

double Geo_Constraint_Dihedral(unsigned n, const double *x, double *grad, void *data)
{
	int ia, ib, ic, id, iPos;
	double phi, phi0, d_phi, g[4][3];
	CMol* pMol;
	GEO_FIX_DIHEDRAL *Geo_Fix_phi;

	Geo_Fix_phi = (GEO_FIX_DIHEDRAL*)data;
	pMol = (CMol*)(Geo_Fix_phi->pMol);
	ia = Geo_Fix_phi->ia;
	ib = Geo_Fix_phi->ib;
	ic = Geo_Fix_phi->ic;
	id = Geo_Fix_phi->id;
	phi0 = Geo_Fix_phi->phi0;

	if(grad)	{
		phi = pMol->Query_Dihedral(ia, ib, ic, id, 1, g);
		memset(grad, 0, sizeof(double)*n);

		iPos = 3 * ia;
		grad[iPos  ] = g[0][0];
		grad[iPos+1] = g[0][1];
		grad[iPos+2] = g[0][2];

		iPos = 3 * ib;
		grad[iPos  ] = g[1][0];
		grad[iPos+1] = g[1][1];
		grad[iPos+2] = g[1][2];

		iPos = 3 * ic;
		grad[iPos  ] = g[2][0];
		grad[iPos+1] = g[2][1];
		grad[iPos+2] = g[2][2];

		iPos = 3 * id;
		grad[iPos  ] = g[3][0];
		grad[iPos+1] = g[3][1];
		grad[iPos+2] = g[3][2];
	}
	else	{
		phi = pMol->Query_Dihedral(ia, ib, ic, id, 0, g);
	}

	d_phi = phi - phi0;
	if(d_phi < -PI)	{
		d_phi += PI2;
	}
	else if(d_phi > PI)	{
		d_phi -= PI2;
	}

	return d_phi;
}

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;
}

double Geo_Constraint_Angle(unsigned n, const double *x, double *grad, void *data)
{
	int ia, ib, ic, iPos;
	double theta, theta0, g[3][3];
	CMol* pMol;
	GEO_FIX_ANGLE *Geo_Fix_theta;

	Geo_Fix_theta = (GEO_FIX_ANGLE*)data;
	pMol = (CMol*)(Geo_Fix_theta->pMol);
	ia = Geo_Fix_theta->ia;
	ib = Geo_Fix_theta->ib;
	ic = Geo_Fix_theta->ic;
	theta0 = Geo_Fix_theta->theta0;

	if(grad)	{
		theta = pMol->Query_Angle(ia, ib, ic, 1, g);
		memset(grad, 0, sizeof(double)*n);

		iPos = 3 * ia;
		grad[iPos  ] = g[0][0];
		grad[iPos+1] = g[0][1];
		grad[iPos+2] = g[0][2];

		iPos = 3 * ib;
		grad[iPos  ] = g[1][0];
		grad[iPos+1] = g[1][1];
		grad[iPos+2] = g[1][2];

		iPos = 3 * ic;
		grad[iPos  ] = g[2][0];
		grad[iPos+1] = g[2][1];
		grad[iPos+2] = g[2][2];
	}
	else	{
		theta = pMol->Query_Angle(ia, ib, ic, 0, g);
	}

	return (theta - theta0);
}

void Get_Representative_14_Atoms(int iDih)
{
	int ia, ib, ic, id, i, Idx, nAtom_Connected, nAtom_Connected_Max=0;

	ib = Dih_Bond[iDih][1];
	ic = Dih_Bond[iDih][2];

	//start	to find the representative atom for the left side, ia
	nAtom_Connected_Max=-100;
	ia = -1;
	for(i=0; i<Bond_Count[ib]; i++)	{
		Idx = Bond_List[ib][i];
		if(Idx == ic)	{
			continue;
		}
		nAtom_Connected = Mol.Count_All_Atoms_Connected(Idx, ib);
		if(nAtom_Connected > nAtom_Connected_Max)	{
			nAtom_Connected_Max = nAtom_Connected;
			ia = Idx;
		}
		else if( (i!=0) && (nAtom_Connected == nAtom_Connected_Max) && (Mol.mass[Idx] > Mol.mass[ia]) )	{
			ia = Idx;
		}
	}
	if(ia>=0)	{
		Dih_Bond[iDih][0] = ia;
	}
	else	{
		Quit_With_Error_Msg("Fail to find the representative atom (heavy atom), ia, in dihedral.\nQuit\n");
	}
	//end	to find the representative atom for the left side, ia


	//start	to find the representative atom for the left side, ia
	nAtom_Connected_Max=-100;
	id = -1;
	for(i=0; i<Bond_Count[ic]; i++)	{
		Idx = Bond_List[ic][i];
		if(Idx == ib)	{
			continue;
		}
		nAtom_Connected = Mol.Count_All_Atoms_Connected(Idx, ic);
		if(nAtom_Connected > nAtom_Connected_Max)	{
			nAtom_Connected_Max = nAtom_Connected;
			id = Idx;
		}
		else if( (i!=0) && (nAtom_Connected == nAtom_Connected_Max) && (Mol.mass[Idx] > Mol.mass[id]) )	{
			id = Idx;
		}
	}
	if(id>=0)	{
		Dih_Bond[iDih][3] = id;
	}
	else	{
		Quit_With_Error_Msg("Fail to find the representative atom (heavy atom), id, in dihedral.\nQuit\n");
	}
	//end	to find the representative atom for the left side, ia

	return;
}

void SaveOptPdb(char szName[])
{
	int i;
	for(i=0; i<n_Phi; i++)	{
		Mol.QueryDihedral(IdxDihSelect[i]);
		Mol.Edit_Dihedral(IdxDihSelect[i], Phi_To_Set_Scan[i]);
	}
	Mol.SavePdb(szName);
}

void To_Setup_All_Dihedral_Constraint(void)
{
	int i;

	Mol_ESP.nPhi_Fixed = 0;

	for(i=0; i<n_Phi; i++)	{
		Mol_ESP.To_Fix_Dihedral(DihList[i][0], DihList[i][1], DihList[i][2], DihList[i][3], 
			Mol_ESP.Query_Dihedral(DihList[i][0], DihList[i][1], DihList[i][2], DihList[i][3], 0, NULL));
	}
}

void Gen_Soft_Dihedral_List_MD_High_T(void)
{
	int Step, Idx, RecIdx, iPos, ia, ib, ic, id;

	Mol.Init_LangevinDynamics(T_Sim);

	nDihedral = Mol.nDihedral;

	for(Step=1; Step<=MAX_STEP; Step++)	{
		Mol.LangevinDynamics(Step);
		if(Step%GAP==0)	{
			RecIdx = Step/GAP - 1;
			for(Idx=0; Idx<nDihedral; Idx++)	{
				dih_List[RecIdx][Idx] = Mol.dih_Phi_List[Idx];
			}
		}
	}

	Cal_Sig();
	Cal_Sig_Shift_PI2();

	n_Soft_Dih = 0;
	for(Idx=0; Idx<nDihedral; Idx++)	{
		if(IsRigidDih[Idx] == 0)	{
			iPos = Idx*4;
			ia = Mol.DihedralList[iPos  ];
			ib = Mol.DihedralList[iPos+1];
			ic = Mol.DihedralList[iPos+2];
			id = Mol.DihedralList[iPos+3];
			if(ib <= ic)	{
				if(!Is_In_Soft_Dihedral_List(ib, ic))	{
					Dih_Bond[n_Soft_Dih][0] = ia;
					Dih_Bond[n_Soft_Dih][1] = ib;
					Dih_Bond[n_Soft_Dih][2] = ic;
					Dih_Bond[n_Soft_Dih][3] = id;
					n_Soft_Dih++;
					printf("%3d %3d %3d %3d \n", Mol.DihedralList[iPos]+1, Mol.DihedralList[iPos+1]+1, Mol.DihedralList[iPos+2]+1, Mol.DihedralList[iPos+3]+1);
				}
			}
			else	{
				if(!Is_In_Soft_Dihedral_List(ic, ib))	{
					Dih_Bond[n_Soft_Dih][0] = id;
					Dih_Bond[n_Soft_Dih][1] = ic;
					Dih_Bond[n_Soft_Dih][2] = ib;
					Dih_Bond[n_Soft_Dih][3] = ia;
					n_Soft_Dih++;
					printf("%3d %3d %3d %3d \n", Mol.DihedralList[iPos]+1, Mol.DihedralList[iPos+1]+1, Mol.DihedralList[iPos+2]+1, Mol.DihedralList[iPos+3]+1);
				}
			}
		}
	}
}

int Read_Soft_DihedralList(void)
{
	FILE *fIn;
	int ReadItem;
	char szLine[256], *ReadLine;

	n_Phi = 0;

	fIn = fopen(szPhiToScan, "r");

	while(1)	{
		if(feof(fIn))	{
			break;
		}
		ReadLine = fgets(szLine, 128, fIn);
		if(ReadLine == NULL)	{
			break;
		}

		ReadItem = sscanf(szLine, "%d %d %d %d %lf %lf %lf %lf %lf %lf", 
			&(DihList[n_Phi][0]), &(DihList[n_Phi][1]), &(DihList[n_Phi][2]), &(DihList[n_Phi][3]), 
			&(Phi_Set[n_Phi][0]), &(Phi_Set[n_Phi][1]), &(Phi_Set[n_Phi][2]), &(Phi_Set[n_Phi][3]), &(Phi_Set[n_Phi][4]), &(Phi_Set[n_Phi][5]));

		if(ReadItem > 4)	{
			DihList[n_Phi][0]--;
			DihList[n_Phi][1]--;
			DihList[n_Phi][2]--;
			DihList[n_Phi][3]--;
			IdxDihSelect[n_Phi] = Mol.Query_Dihedral_Index(DihList[n_Phi][0], DihList[n_Phi][1], DihList[n_Phi][2], DihList[n_Phi][3]);

			if(IdxDihSelect[n_Phi] < 0)	{
				Quit_With_Error_Msg("Fail to identify the index of one soft dihedral.\n");
			}
			Mol.BuildSegmentList_Dihedrals(IdxDihSelect[n_Phi]);

			n_State_List[n_Phi] = ReadItem - 4;
			n_Phi++;
		}
		else	{
			break;
		}
	}

	fclose(fIn);

	memcpy(Phi_Set_Save, Phi_Set, sizeof(double)*N_MAX_DIH*MAX_N_STATE);
	memcpy(n_State_List_Save, n_State_List, sizeof(int)*N_MAX_DIH);

	return n_Phi;
}

void Output_Gaussian_File(void)
{
	FILE *fOut, *fIn;
	int i, nAtom, IdxScan;
	char szNameGjf[256], szCmd[256], szName_Output[256];

	Determine_QM_Scan_Part();

	nAtom = Mol.nAtom;
	for(IdxScan=0; IdxScan<n_QM_Part; IdxScan++)	{
		Mol.QueryDihedral(IdxDihSelect[Active_Phi]);
		Mol.Edit_Dihedral(IdxDihSelect[Active_Phi], Phi_QM_Start[IdxScan]);

		sprintf(szNameGjf, "qm-scan-%d-%d.gjf", Active_Phi+1, IdxScan+1);
		fOut = fopen(szNameGjf, "w");
		fprintf(fOut, "%s", szQM_Level_1D_Scan_Cmd);
		for(i=0; i<nAtom; i++)	{
//			fprintf(fOut, "%c  %9.5lf %9.5lf %9.5lf\n", Mol.AtomName[i][0], Mol.x[i], Mol.y[i], Mol.z[i]);
			fprintf(fOut, "%s  %9.5lf %9.5lf %9.5lf\n", szElemName[i], Mol.x[i], Mol.y[i], Mol.z[i]);
		}
		fprintf(fOut, "\n");
		
		for(i=0; i<n_Phi; i++)	{
			if(i == Active_Phi)	{
				fprintf(fOut, "%s", szQM_Scan[IdxScan]);
			}
			else	{
				fprintf(fOut, "%d %d %d %d F\n", DihList[i][0]+1, DihList[i][1]+1, DihList[i][2]+1, DihList[i][3]+1);
			}
		}
		fprintf(fOut, "\n");
		fclose(fOut);

		sprintf(szName_Output, "qm-1d-phi-%d-p%d.out", Active_Phi+1, IdxScan+1);
		sprintf(szCmd, "%s < %s > %s", szExe_G09, szNameGjf, szName_Output);


		fIn = fopen(szName_Output, "r");
		if(fIn == NULL)	{
			system(szCmd);
		}
		else	{
			fclose(fIn);
		}

		Extract_Coord_E(szName_Output, Active_Phi, IdxScan);
	}
}

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;
}

int main(int argc, char *argv[])
{
	if(argc < 4)	{
		printf("Usage: add-tip3 your-rtf your-prm your-psf E14FAC\nQuit\n");
		return 1;
	}

	strcpy(szRtfFile, argv[1]);
	strcpy(szPrmFile, argv[2]);
	strcpy(szXpsfName, argv[3]);

	if(argc == 5)	{
		E14FAC = atof(argv[4]);
	}

	fFile_Run_Log = fopen("log-modify-rtf.txt", "w");
	Mol.ReadPSF(szXpsfName, 0);
	GetAtomTypeFromMass();
	Count_Bond_and_H();
	Assign_Alpha_Miller();

	Add_Tip3_Rtf_File();

	Add_Tip3_Prm_File();

	fclose(fFile_Run_Log);

	return 0;
}

int main(int argc, char *argv[])
{
	if(argc == 2)	{
		To_Exclude_Methyl_Torsion = 0;
	}
	else	{
		To_Exclude_Methyl_Torsion = 1;
	}

	memset(IsRigidDih, 0, sizeof(int)*MAX_DIHEDRAL_LIST);

	fFile_Run_Log = fopen("identify-soft.log", "w");

	ForceField.ReadForceField(szForceFiled);
	Mol.ReadPSF(szXpsfFile, 0);
	Mol.AssignForceFieldParameters(&ForceField);
	Mol.ReadCRD(szCrdFile);
	Setup_Bond_List();

//	Gen_Soft_Dihedral_List_MD_High_T();
//	printf("\n\n");

	Gen_Soft_Dihedral_List_Excluding_Ring();
	printf("\n\n");

	To_Refine_Soft_Dihedral_List();

/*	
	fData = fopen("dihedral.txt", "w");
	for(Step=0; Step<MAX_STEP/GAP; Step++)	{
		fprintf(fData, "%8d ", Step);
		for(Idx=0; Idx<nDihedral; Idx++)	{
			if(IsRigidDih[Idx] == 0)	{
				fprintf(fData, " %7.2lf", dih_List[Step][Idx]);
			}
		}
		fprintf(fData, "\n");
	}
	fclose(fData);
*/
	fclose(fFile_Run_Log);


	return 0;
}

int Read_Soft_DihedralList(void)
{
    FILE *fIn;
    int ReadItem;
    char szLine[256], *ReadLine, ErrorMsg[256], szFlag[256];

    n_Phi = 0;

    fIn = fopen(szPhiToScan, "r");

    while(1)	{
        if(feof(fIn))	{
            sprintf(ErrorMsg, "Fail to open file: %s\nQuit\n", szPhiToScan);
            Quit_With_Error_Msg(ErrorMsg);
        }
        ReadLine = fgets(szLine, 128, fIn);
        if(ReadLine == NULL)	{
            break;
        }

        ReadItem = sscanf(szLine, "%d %d %d %d %s",
                          &(DihList[n_Phi][0]), &(DihList[n_Phi][1]), &(DihList[n_Phi][2]), &(DihList[n_Phi][3]), szFlag);

        if(ReadItem >= 4)	{
            DihList[n_Phi][0]--;
            DihList[n_Phi][1]--;
            DihList[n_Phi][2]--;
            DihList[n_Phi][3]--;
            IdxDihSelect[n_Phi] = Mol_ESP.Query_Dihedral_Index(DihList[n_Phi][0], DihList[n_Phi][1], DihList[n_Phi][2], DihList[n_Phi][3]);

            if(IdxDihSelect[n_Phi] < 0)	{
                Quit_With_Error_Msg("Fail to identify the index of one soft dihedral.\n");
            }

            Is_Dih_Fixed[n_Phi] = 0;
            if(ReadItem==5)	{
                if( (strcmp(szFlag, "F")==0) || (strcmp(szFlag, "f")==0) )	{
                    Is_Dih_Fixed[n_Phi] = 1;
                }
            }

            n_Phi++;
        }
        else	{
            break;
        }
    }

    fclose(fIn);

    Setup_Phi_Constraint_List();

    return n_Phi;
}

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;
}