Skip to content

Commit

Permalink
Merge pull request #15 from jnochebuena/main
Browse files Browse the repository at this point in the history 
Geometry optimization using GEM
  • Loading branch information
@g-andres
g-andres authored Aug 8, 2024
2 parents d378f2c + 879ddfc commit 2ba4d38
Show file tree
Hide file tree
Showing 2 changed files with 46 additions and 18 deletions.
6 changes: 4 additions & 2 deletions src/Struct_writer.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -751,7 +751,8 @@ void WritePSITHONInput(vector<QMMMAtom>& QMMMData, string calcTyp,
call << '"' << '"' << '"' << ")" << '\n';

call << "molB = psi4.geometry(" << '"' << '"' << '"' << '\n';
call << " " << QMMMOpts.charge;
//call << " " << QMMMOpts.charge;
call << " " << "0"; //zero by convenience, fix later
call << " " << QMMMOpts.spin << '\n';
for (int i=0;i<Natoms;i++)
{
Expand All @@ -776,7 +777,7 @@ void WritePSITHONInput(vector<QMMMAtom>& QMMMData, string calcTyp,
call << "K_exchange = " << QMMMOpts.kexchange << "\n" << "\n";

call << "epsilon = 1e-12" << '\n';
call << "regularizer = 0.000001" << '\n';
call << "regularizer = 1e-8" << '\n';
call << "constrain = True" << '\n' << '\n';

call << "psi4.set_options({'puream' : False, 'print' : 1,'scf_type' : 'df'})" << '\n' << '\n';
Expand Down Expand Up @@ -855,6 +856,7 @@ void WritePSITHONInput(vector<QMMMAtom>& QMMMData, string calcTyp,
call << " fit_coefficients += lam * qdot" << '\n';
call << " new_nelec = np.dot(fit_coefficients, q)" << '\n';
call << " psi4.core.print_out(f'\\nNumber of electrons after constraining: {new_nelec:.6f}')" << '\n';
call << " np.savetxt('fitcoeffs.dat',fit_coefficients)" << '\n';
call << " else:" << '\n';
call << " fit_coefficients = np.loadtxt('fitcoeffs.dat')" << '\n';
call << " coefs_vec = psi4.core.Vector(aux_basis.nbf())" << '\n';
Expand Down
58 changes: 42 additions & 16 deletions src/TINKER.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -873,10 +873,16 @@ double TINKERForces(vector<QMMMAtom>& QMMMData, VectorXd& forces,
else{
call.str("");
call << "rm -f";
call << " LICHM_" << bead << ".xyz";
call << " LICHM_" << bead << ".key";
//call << " LICHM_" << bead << ".xyz";
//call << " LICHM_" << bead << ".key";
call << " LICHM_" << bead << ".err";
call << "; mv ";
call << "LICHM_" << bead << ".xyz ";
call << "LICHM_TINKERForces_" << bead << ".xyz ";
call << "; mv ";
call << "LICHM_" << bead << ".key ";
call << "LICHM_TINKERForces_" << bead << ".key ";
call << "; mv ";
call << "LICHM_" << bead << ".grad ";
call << "LICHM_TINKERForces_" << bead << ".grad ";
globalSys = system(call.str().c_str());
Expand Down Expand Up @@ -1429,10 +1435,16 @@ double TINKERPolForces(vector<QMMMAtom>& QMMMData, VectorXd& forces,
else{
call.str("");
call << "rm -f";
call << " LICHM_" << bead << ".xyz";
call << " LICHM_" << bead << ".key";
//call << " LICHM_" << bead << ".xyz";
//call << " LICHM_" << bead << ".key";
call << " LICHM_" << bead << ".err";
call << "; mv ";
call << "LICHM_" << bead << ".xyz ";
call << "LICHM_TINKERPolForces_" << bead << ".xyz ";
call << "; mv ";
call << "LICHM_" << bead << ".key ";
call << "LICHM_TINKERPolForces_" << bead << ".key ";
call << "; mv ";
call << "LICHM_" << bead << ".grad ";
call << "LICHM_TINKERPolForces_" << bead << ".grad ";
globalSys = system(call.str().c_str());
Expand Down Expand Up @@ -1481,6 +1493,7 @@ double TINKEREnergy(vector<QMMMAtom>& QMMMData, QMMMSettings& QMMMOpts,
{
outFile << "ONLY-HALGREN" << '\n';
}
//E:JORGE
if (QMMMOpts.useLREC)
{
//Apply cutoff
Expand Down Expand Up @@ -1642,9 +1655,6 @@ double TINKEREnergy(vector<QMMMAtom>& QMMMData, QMMMSettings& QMMMOpts,
inFile.open(call.str().c_str(),ios_base::in);
//Read MM potential energy
bool EFound = 0;
//S:JORGE
//cout << "==> " << E << endl;
//E:JORGE
while ((!inFile.eof()) and inFile.good())
{
inFile >> dummy;
Expand All @@ -1658,9 +1668,6 @@ double TINKEREnergy(vector<QMMMAtom>& QMMMData, QMMMSettings& QMMMOpts,
}
}
}
//S:JORGE
//cout << "==> " << E << endl;
//E:JORGE
if (!EFound)
{
//Warn user if no energy was found
Expand Down Expand Up @@ -2234,6 +2241,12 @@ double TINKEROpt(vector<QMMMAtom>& QMMMData, QMMMSettings& QMMMOpts, int bead,
outFile << "#LICHEM MM keywords"; //Marks the changes
}
outFile << '\n';
//S:JORGE
if (GEM)
{
outFile << "ONLY-HALGREN" << '\n';
}
//E:JORGE
if (QMMMOpts.useMMCut)
{
//Apply cutoff
Expand Down Expand Up @@ -2373,7 +2386,7 @@ double TINKEROpt(vector<QMMMAtom>& QMMMData, QMMMSettings& QMMMOpts, int bead,
}
}
}
if (AMOEBA)
if (AMOEBA or GEM)
{
for (int i=0;i<Natoms;i++)
{
Expand Down Expand Up @@ -2460,6 +2473,13 @@ double TINKEROpt(vector<QMMMAtom>& QMMMData, QMMMSettings& QMMMOpts, int bead,
}
outFile.flush();
outFile.close();
//S:JORGE
//Remove previous structures before minimizing
call.str("");
call << "rm -f";
call << " LICHM_" << bead << ".xyz_*";
globalSys = system(call.str().c_str());
//E:JORGE
//Run optimization
call.str("");
call << "minimize LICHM_";
Expand All @@ -2472,7 +2492,7 @@ double TINKEROpt(vector<QMMMAtom>& QMMMData, QMMMSettings& QMMMOpts, int bead,
call << "LICHM_" << bead << ".xyz_2";
inFile.open(call.str().c_str(),ios_base::in);
getline(inFile,dummy); //Discard number of atoms
if (PBCon)
if (PBCon and !GEM) //Modified for Tinker-HP compatibility
{
//Discard PBC information
getline(inFile,dummy);
Expand All @@ -2494,19 +2514,25 @@ double TINKEROpt(vector<QMMMAtom>& QMMMData, QMMMSettings& QMMMOpts, int bead,
call << "rm -f";
call << " LICHM_" << bead << ".xyz";
call << " LICHM_" << bead << ".log";
call << " LICHM_" << bead << ".xyz_*";
//call << " LICHM_" << bead << ".xyz_*";
call << " LICHM_" << bead << ".key";
call << " LICHM_" << bead << ".err";
globalSys = system(call.str().c_str());
}
else{
call.str("");
call << "rm -f";
call << " LICHM_" << bead << ".xyz";
call << " LICHM_" << bead << ".xyz_*";
call << " LICHM_" << bead << ".key";
//call << " LICHM_" << bead << ".xyz";
//call << " LICHM_" << bead << ".xyz_*";
//call << " LICHM_" << bead << ".key";
call << " LICHM_" << bead << ".err";
call << "; mv ";
call << "LICHM_" << bead << ".xyz ";
call << "LICHM_TINKEROpt_" << bead << ".xyz ";
call << "; mv ";
call << "LICHM_" << bead << ".key ";
call << "LICHM_TINKEROpt_" << bead << ".key ";
call << "; mv ";
call << "LICHM_" << bead << ".log ";
call << "LICHM_TINKEROpt_" << bead << ".log ";
globalSys = system(call.str().c_str());
Expand Down

0 comments on commit 2ba4d38

Please sign in to comment.