From 1ffb0a711232f353283d5043eb65b1817c040031 Mon Sep 17 00:00:00 2001 From: mohammad-ghahremanpour Date: Sat, 11 Jan 2020 22:58:37 -0500 Subject: [PATCH] changes for fitting atmomic polarizability in tune_eem Change-Id: I90fe5c64659dd060afdb27bb745aa94efc21926e --- src/programs/alexandria/mymol.cpp | 483 +++++++++++++++++----------------- src/programs/alexandria/mymol_low.cpp | 20 ++ src/programs/alexandria/mymol_low.h | 8 + src/programs/alexandria/poldata.cpp | 97 +++++++ src/programs/alexandria/poldata.h | 21 +- src/programs/alexandria/tune_eem.cpp | 18 +- src/programs/alexandria/tune_zeta.cpp | 6 +- 7 files changed, 400 insertions(+), 253 deletions(-) diff --git a/src/programs/alexandria/mymol.cpp b/src/programs/alexandria/mymol.cpp index 1cf85f4b4e..ec05f23db9 100644 --- a/src/programs/alexandria/mymol.cpp +++ b/src/programs/alexandria/mymol.cpp @@ -2200,295 +2200,288 @@ void MyMol::UpdateIdef(const Poldata *pd, nonbondedFromPdToMtop(mtop_, atoms_, pd, fr_); if (debug) { - pr_ffparams(debug, 0, "UpdateIdef Before", &mtop_->ffparams, FALSE); + pr_ffparams(debug, 0, "UpdateIdef Before", &mtop_->ffparams, false); } - return; } - auto fs = pd->findForces(iType); - if (fs == pd->forcesEnd()) + else if (iType == eitPOLARIZATION) { - gmx_fatal(FARGS, "Can not find the force %s to update", - iType2string(iType)); - } - auto ftype = fs->fType(); - // Make a list of bonded types that can be indexed - // with the atomtype. That should speed up the code - // below somewhat. - std::vector btype(mtop_->ffparams.atnr); - for (int i = 0; i < mtop_->natoms; i++) - { - std::string bt; - if (pd->atypeToBtype(*atoms_->atomtype[i], bt)) + auto pw = SearchPlist(plist_, iType); + if (plist_.end() != pw) { - btype[atoms_->atom[i].type] = bt; - } - else if (!(atoms_->atom[i].ptype == eptShell || - atoms_->atom[i].ptype == eptVSite)) - { - gmx_fatal(FARGS, "Cannot find bonded type for atomtype %s", - *atoms_->atomtype[i]); + auto ft = pw->getFtype(); + polarizabilityFromPdToMtop(mtop_, ltop_, atoms_, pd, ft); } } - switch (iType) + else { - case eitBONDS: + // Update other iTypes + auto fs = pd->findForces(iType); + if (fs == pd->forcesEnd()) + { + gmx_fatal(FARGS, "Can not find the force %s to update", + iType2string(iType)); + } + auto ftype = fs->fType(); + // Make a list of bonded types that can be indexed + // with the atomtype. That should speed up the code + // below somewhat. + std::vector btype(mtop_->ffparams.atnr); + for (int i = 0; i < mtop_->natoms; i++) { - auto lu = string2unit(fs->unit().c_str()); - if (-1 == lu) + std::string bt; + if (pd->atypeToBtype(*atoms_->atomtype[i], bt)) { - gmx_fatal(FARGS, "Unknown length unit '%s' for bonds", - fs->unit().c_str()); + btype[atoms_->atom[i].type] = bt; } - for (auto i = 0; i < ltop_->idef.il[ftype].nr; i += interaction_function[ftype].nratoms+1) + else if (!(atoms_->atom[i].ptype == eptShell || + atoms_->atom[i].ptype == eptVSite)) { - auto tp = ltop_->idef.il[ftype].iatoms[i]; - std::vector atoms; - for (int j = 1; j < interaction_function[ftype].nratoms+1; j++) - { - atoms.push_back(btype[atoms_->atom[ltop_->idef.il[ftype].iatoms[i+j]].type]); - } - double value, sigma; - size_t ntrain; - if (pd->searchForceIType(atoms, params, &value, &sigma, &ntrain, iType)) - { - auto bondLength = convert2gmx(value, lu); - auto parameters = gmx::splitString(params); - switch (ftype) - { - case F_MORSE: - { - mtop_->ffparams.iparams[tp].morse.b0A = - mtop_->ffparams.iparams[tp].morse.b0B = - ltop_->idef.iparams[tp].morse.b0A = - ltop_->idef.iparams[tp].morse.b0B = bondLength; - GMX_RELEASE_ASSERT(parameters.size() == 2, "Need exactly two parameters for Morse bonds"); - mtop_->ffparams.iparams[tp].morse.cbA = - mtop_->ffparams.iparams[tp].morse.cbB = - ltop_->idef.iparams[tp].morse.cbA = - ltop_->idef.iparams[tp].morse.cbB = - gmx::doubleFromString(parameters[0].c_str()); - mtop_->ffparams.iparams[tp].morse.betaA = - mtop_->ffparams.iparams[tp].morse.betaB = - ltop_->idef.iparams[tp].morse.betaA = - ltop_->idef.iparams[tp].morse.betaB = - gmx::doubleFromString(parameters[1].c_str()); - } - break; - default: - gmx_fatal(FARGS, "Unsupported bondtype %s in UpdateIdef", - fs->function().c_str()); - } - } + gmx_fatal(FARGS, "Cannot find bonded type for atomtype %s", + *atoms_->atomtype[i]); } } - break; - case eitANGLES: - case eitLINEAR_ANGLES: + switch (iType) { - for (auto i = 0; i < ltop_->idef.il[ftype].nr; i += interaction_function[ftype].nratoms+1) + case eitBONDS: { - auto tp = ltop_->idef.il[ftype].iatoms[i]; - std::vector atoms; - for (int j = 1; j < interaction_function[ftype].nratoms+1; j++) + auto lu = string2unit(fs->unit().c_str()); + if (-1 == lu) { - atoms.push_back(btype[atoms_->atom[ltop_->idef.il[ftype].iatoms[i+j]].type]); + gmx_fatal(FARGS, "Unknown length unit '%s' for bonds", + fs->unit().c_str()); } - double angle, sigma; - size_t ntrain; - if (pd->searchForceIType(atoms, params, &angle, &sigma, &ntrain, iType)) + for (auto i = 0; i < ltop_->idef.il[ftype].nr; i += interaction_function[ftype].nratoms+1) { - auto parameters = gmx::splitString(params); - auto r13 = calc_r13(pd, atoms[0], atoms[1], atoms[2], angle); - switch (ftype) + auto tp = ltop_->idef.il[ftype].iatoms[i]; + std::vector atoms; + for (int j = 1; j < interaction_function[ftype].nratoms+1; j++) { - case F_ANGLES: - { - mtop_->ffparams.iparams[tp].harmonic.rA = - mtop_->ffparams.iparams[tp].harmonic.rB = - ltop_->idef.iparams[tp].harmonic.rA = - ltop_->idef.iparams[tp].harmonic.rB = angle; - GMX_RELEASE_ASSERT(parameters.size() == 1, "Need exactly one parameters for Harmonic angles"); - mtop_->ffparams.iparams[tp].harmonic.krA = - mtop_->ffparams.iparams[tp].harmonic.krB = - ltop_->idef.iparams[tp].harmonic.krA = - ltop_->idef.iparams[tp].harmonic.krB = - gmx::doubleFromString(parameters[0].c_str()); - } - break; - case F_UREY_BRADLEY: - { - mtop_->ffparams.iparams[tp].u_b.thetaA = - mtop_->ffparams.iparams[tp].u_b.thetaB = - ltop_->idef.iparams[tp].u_b.thetaA = - ltop_->idef.iparams[tp].u_b.thetaB = angle; - - mtop_->ffparams.iparams[tp].u_b.r13A = - mtop_->ffparams.iparams[tp].u_b.r13B = - ltop_->idef.iparams[tp].u_b.r13A = - ltop_->idef.iparams[tp].u_b.r13B = r13; - GMX_RELEASE_ASSERT(parameters.size() == 2, "Need exactly two parameters for Urey-Bradley angles"); - - mtop_->ffparams.iparams[tp].u_b.kthetaA = - mtop_->ffparams.iparams[tp].u_b.kthetaB = - ltop_->idef.iparams[tp].u_b.kthetaA = - ltop_->idef.iparams[tp].u_b.kthetaB = - gmx::doubleFromString(parameters[0].c_str()); - mtop_->ffparams.iparams[tp].u_b.kUBA = - mtop_->ffparams.iparams[tp].u_b.kUBB = - ltop_->idef.iparams[tp].u_b.kUBA = - ltop_->idef.iparams[tp].u_b.kUBB = - gmx::doubleFromString(parameters[1].c_str()); - } - break; - case F_LINEAR_ANGLES: + atoms.push_back(btype[atoms_->atom[ltop_->idef.il[ftype].iatoms[i+j]].type]); + } + double value, sigma; + size_t ntrain; + if (pd->searchForceIType(atoms, params, &value, &sigma, &ntrain, iType)) + { + auto bondLength = convert2gmx(value, lu); + auto parameters = gmx::splitString(params); + switch (ftype) { - double relative_position = calc_relposition(pd, atoms[0], atoms[1], atoms[2]); - - mtop_->ffparams.iparams[tp].linangle.aA = - mtop_->ffparams.iparams[tp].linangle.aB = - ltop_->idef.iparams[tp].linangle.aA = - ltop_->idef.iparams[tp].linangle.aB = relative_position; - - mtop_->ffparams.iparams[tp].linangle.r13A = - mtop_->ffparams.iparams[tp].linangle.r13B = - ltop_->idef.iparams[tp].linangle.r13A = - ltop_->idef.iparams[tp].linangle.r13B = r13; - GMX_RELEASE_ASSERT(parameters.size() == 2, "Need exactly two parameters for Linear angles"); - - mtop_->ffparams.iparams[tp].linangle.klinA = - mtop_->ffparams.iparams[tp].linangle.klinB = - ltop_->idef.iparams[tp].linangle.klinA = - ltop_->idef.iparams[tp].linangle.klinB = - gmx::doubleFromString(parameters[0].c_str()); - mtop_->ffparams.iparams[tp].linangle.kUBA = - mtop_->ffparams.iparams[tp].linangle.kUBB = - ltop_->idef.iparams[tp].linangle.kUBA = - ltop_->idef.iparams[tp].linangle.kUBB = - gmx::doubleFromString(parameters[1].c_str()); + case F_MORSE: + { + mtop_->ffparams.iparams[tp].morse.b0A = + mtop_->ffparams.iparams[tp].morse.b0B = + ltop_->idef.iparams[tp].morse.b0A = + ltop_->idef.iparams[tp].morse.b0B = bondLength; + GMX_RELEASE_ASSERT(parameters.size() == 2, "Need exactly two parameters for Morse bonds"); + mtop_->ffparams.iparams[tp].morse.cbA = + mtop_->ffparams.iparams[tp].morse.cbB = + ltop_->idef.iparams[tp].morse.cbA = + ltop_->idef.iparams[tp].morse.cbB = + gmx::doubleFromString(parameters[0].c_str()); + mtop_->ffparams.iparams[tp].morse.betaA = + mtop_->ffparams.iparams[tp].morse.betaB = + ltop_->idef.iparams[tp].morse.betaA = + ltop_->idef.iparams[tp].morse.betaB = + gmx::doubleFromString(parameters[1].c_str()); + } + break; + default: + gmx_fatal(FARGS, "Unsupported bondtype %s in UpdateIdef", + fs->function().c_str()); } - break; - default: - gmx_fatal(FARGS, "Unsupported angletype %s in UpdateIdef", - fs->function().c_str()); } } } - } - break; - case eitPROPER_DIHEDRALS: - case eitIMPROPER_DIHEDRALS: - { - for (auto i = 0; i < ltop_->idef.il[ftype].nr; i += interaction_function[ftype].nratoms+1) + break; + case eitANGLES: + case eitLINEAR_ANGLES: { - auto tp = ltop_->idef.il[ftype].iatoms[i]; - std::vector atoms; - for (int j = 1; j < interaction_function[ftype].nratoms+1; j++) - { - atoms.push_back(btype[atoms_->atom[ltop_->idef.il[ftype].iatoms[i+j]].type]); - } - double angle, sigma; - size_t ntrain; - if (pd->searchForceIType(atoms, params, - &angle, &sigma, &ntrain, iType)) + for (auto i = 0; i < ltop_->idef.il[ftype].nr; i += interaction_function[ftype].nratoms+1) { - auto parameters = gmx::splitString(params); - switch (ftype) + auto tp = ltop_->idef.il[ftype].iatoms[i]; + std::vector atoms; + for (int j = 1; j < interaction_function[ftype].nratoms+1; j++) { - case F_FOURDIHS: + atoms.push_back(btype[atoms_->atom[ltop_->idef.il[ftype].iatoms[i+j]].type]); + } + double angle, sigma; + size_t ntrain; + if (pd->searchForceIType(atoms, params, &angle, &sigma, &ntrain, iType)) + { + auto parameters = gmx::splitString(params); + auto r13 = calc_r13(pd, atoms[0], atoms[1], atoms[2], angle); + switch (ftype) { - std::vector old; - for (auto parm : parameters) + case F_ANGLES: { - old.push_back(gmx::doubleFromString(parm.c_str())); + mtop_->ffparams.iparams[tp].harmonic.rA = + mtop_->ffparams.iparams[tp].harmonic.rB = + ltop_->idef.iparams[tp].harmonic.rA = + ltop_->idef.iparams[tp].harmonic.rB = angle; + GMX_RELEASE_ASSERT(parameters.size() == 1, "Need exactly one parameters for Harmonic angles"); + mtop_->ffparams.iparams[tp].harmonic.krA = + mtop_->ffparams.iparams[tp].harmonic.krB = + ltop_->idef.iparams[tp].harmonic.krA = + ltop_->idef.iparams[tp].harmonic.krB = + gmx::doubleFromString(parameters[0].c_str()); } - auto newparam = &mtop_->ffparams.iparams[tp]; - newparam->rbdihs.rbcA[0] = old[1]+0.5*(old[0]+old[2]); - newparam->rbdihs.rbcA[1] = 0.5*(3.0*old[2]-old[0]); - newparam->rbdihs.rbcA[2] = 4.0*old[3]-old[1]; - newparam->rbdihs.rbcA[3] = -2.0*old[2]; - newparam->rbdihs.rbcA[4] = -4.0*old[3]; - newparam->rbdihs.rbcA[5] = 0.0; - for (int k = 0; k < NR_RBDIHS; k++) + break; + case F_UREY_BRADLEY: { - ltop_->idef.iparams[tp].rbdihs.rbcA[k] = - ltop_->idef.iparams[tp].rbdihs.rbcB[k] = - newparam->rbdihs.rbcB[k] = - newparam->rbdihs.rbcA[k]; + mtop_->ffparams.iparams[tp].u_b.thetaA = + mtop_->ffparams.iparams[tp].u_b.thetaB = + ltop_->idef.iparams[tp].u_b.thetaA = + ltop_->idef.iparams[tp].u_b.thetaB = angle; + + mtop_->ffparams.iparams[tp].u_b.r13A = + mtop_->ffparams.iparams[tp].u_b.r13B = + ltop_->idef.iparams[tp].u_b.r13A = + ltop_->idef.iparams[tp].u_b.r13B = r13; + GMX_RELEASE_ASSERT(parameters.size() == 2, "Need exactly two parameters for Urey-Bradley angles"); + + mtop_->ffparams.iparams[tp].u_b.kthetaA = + mtop_->ffparams.iparams[tp].u_b.kthetaB = + ltop_->idef.iparams[tp].u_b.kthetaA = + ltop_->idef.iparams[tp].u_b.kthetaB = + gmx::doubleFromString(parameters[0].c_str()); + mtop_->ffparams.iparams[tp].u_b.kUBA = + mtop_->ffparams.iparams[tp].u_b.kUBB = + ltop_->idef.iparams[tp].u_b.kUBA = + ltop_->idef.iparams[tp].u_b.kUBB = + gmx::doubleFromString(parameters[1].c_str()); } break; + case F_LINEAR_ANGLES: + { + double relative_position = calc_relposition(pd, atoms[0], atoms[1], atoms[2]); + + mtop_->ffparams.iparams[tp].linangle.aA = + mtop_->ffparams.iparams[tp].linangle.aB = + ltop_->idef.iparams[tp].linangle.aA = + ltop_->idef.iparams[tp].linangle.aB = relative_position; + + mtop_->ffparams.iparams[tp].linangle.r13A = + mtop_->ffparams.iparams[tp].linangle.r13B = + ltop_->idef.iparams[tp].linangle.r13A = + ltop_->idef.iparams[tp].linangle.r13B = r13; + GMX_RELEASE_ASSERT(parameters.size() == 2, "Need exactly two parameters for Linear angles"); + + mtop_->ffparams.iparams[tp].linangle.klinA = + mtop_->ffparams.iparams[tp].linangle.klinB = + ltop_->idef.iparams[tp].linangle.klinA = + ltop_->idef.iparams[tp].linangle.klinB = + gmx::doubleFromString(parameters[0].c_str()); + mtop_->ffparams.iparams[tp].linangle.kUBA = + mtop_->ffparams.iparams[tp].linangle.kUBB = + ltop_->idef.iparams[tp].linangle.kUBA = + ltop_->idef.iparams[tp].linangle.kUBB = + gmx::doubleFromString(parameters[1].c_str()); + } + break; + default: + gmx_fatal(FARGS, "Unsupported angletype %s in UpdateIdef", + fs->function().c_str()); } - case F_PDIHS: - { - mtop_->ffparams.iparams[tp].pdihs.phiA = - mtop_->ffparams.iparams[tp].pdihs.phiB = - ltop_->idef.iparams[tp].pdihs.phiA = - ltop_->idef.iparams[tp].pdihs.phiB = angle; - - GMX_RELEASE_ASSERT(parameters.size() == 2, "Need exactly two parameters for proper dihedrals"); - mtop_->ffparams.iparams[tp].pdihs.cpA = - mtop_->ffparams.iparams[tp].pdihs.cpB = - ltop_->idef.iparams[tp].pdihs.cpA = - ltop_->idef.iparams[tp].pdihs.cpB = - gmx::doubleFromString(parameters[0].c_str()); - /* Multiplicity for Proper Dihedral must be integer - This assumes that the second parameter is Multiplicity */ - mtop_->ffparams.iparams[tp].pdihs.mult = - ltop_->idef.iparams[tp].pdihs.mult = - atoi(parameters[1].c_str()); - } - break; - case F_IDIHS: - { - mtop_->ffparams.iparams[tp].harmonic.rA = - mtop_->ffparams.iparams[tp].harmonic.rB = - ltop_->idef.iparams[tp].harmonic.rA = - ltop_->idef.iparams[tp].harmonic.rB = angle; - - GMX_RELEASE_ASSERT(parameters.size() == 1, "Need exactly one parameter for proper dihedrals"); - mtop_->ffparams.iparams[tp].harmonic.krA = - mtop_->ffparams.iparams[tp].harmonic.krB = - ltop_->idef.iparams[tp].harmonic.krA = - ltop_->idef.iparams[tp].harmonic.krB = - gmx::doubleFromString(parameters[0].c_str()); - } - break; - default: - gmx_fatal(FARGS, "Unsupported dihedral type %s in UpdateIdef", - fs->function().c_str()); } } } - } - break; - case eitLJ14: - case eitPOLARIZATION: - { - auto pw = SearchPlist(plist_, iType); - if (plist_.end() != pw) + break; + case eitPROPER_DIHEDRALS: + case eitIMPROPER_DIHEDRALS: { - auto ft = pw->getFtype(); - for (auto i = 0; i < ltop_->idef.il[ft].nr; i += interaction_function[ft].nratoms+1) + for (auto i = 0; i < ltop_->idef.il[ftype].nr; i += interaction_function[ftype].nratoms+1) { - auto tp = ltop_->idef.il[ft].iatoms[i]; - auto ai = ltop_->idef.il[ft].iatoms[i+1]; - double alpha, sigma; - if (pd->getAtypePol(*atoms_->atomtype[ai], &alpha, &sigma)) + auto tp = ltop_->idef.il[ftype].iatoms[i]; + std::vector atoms; + for (int j = 1; j < interaction_function[ftype].nratoms+1; j++) + { + atoms.push_back(btype[atoms_->atom[ltop_->idef.il[ftype].iatoms[i+j]].type]); + } + double angle, sigma; + size_t ntrain; + if (pd->searchForceIType(atoms, params, + &angle, &sigma, &ntrain, iType)) { - mtop_->ffparams.iparams[tp].polarize.alpha = - ltop_->idef.iparams[tp].polarize.alpha = alpha; + auto parameters = gmx::splitString(params); + switch (ftype) + { + case F_FOURDIHS: + { + std::vector old; + for (auto parm : parameters) + { + old.push_back(gmx::doubleFromString(parm.c_str())); + } + auto newparam = &mtop_->ffparams.iparams[tp]; + newparam->rbdihs.rbcA[0] = old[1]+0.5*(old[0]+old[2]); + newparam->rbdihs.rbcA[1] = 0.5*(3.0*old[2]-old[0]); + newparam->rbdihs.rbcA[2] = 4.0*old[3]-old[1]; + newparam->rbdihs.rbcA[3] = -2.0*old[2]; + newparam->rbdihs.rbcA[4] = -4.0*old[3]; + newparam->rbdihs.rbcA[5] = 0.0; + for (int k = 0; k < NR_RBDIHS; k++) + { + ltop_->idef.iparams[tp].rbdihs.rbcA[k] = + ltop_->idef.iparams[tp].rbdihs.rbcB[k] = + newparam->rbdihs.rbcB[k] = + newparam->rbdihs.rbcA[k]; + } + break; + } + case F_PDIHS: + { + mtop_->ffparams.iparams[tp].pdihs.phiA = + mtop_->ffparams.iparams[tp].pdihs.phiB = + ltop_->idef.iparams[tp].pdihs.phiA = + ltop_->idef.iparams[tp].pdihs.phiB = angle; + + GMX_RELEASE_ASSERT(parameters.size() == 2, "Need exactly two parameters for proper dihedrals"); + mtop_->ffparams.iparams[tp].pdihs.cpA = + mtop_->ffparams.iparams[tp].pdihs.cpB = + ltop_->idef.iparams[tp].pdihs.cpA = + ltop_->idef.iparams[tp].pdihs.cpB = + gmx::doubleFromString(parameters[0].c_str()); + /* Multiplicity for Proper Dihedral must be integer + This assumes that the second parameter is Multiplicity */ + mtop_->ffparams.iparams[tp].pdihs.mult = + ltop_->idef.iparams[tp].pdihs.mult = + atoi(parameters[1].c_str()); + } + break; + case F_IDIHS: + { + mtop_->ffparams.iparams[tp].harmonic.rA = + mtop_->ffparams.iparams[tp].harmonic.rB = + ltop_->idef.iparams[tp].harmonic.rA = + ltop_->idef.iparams[tp].harmonic.rB = angle; + + GMX_RELEASE_ASSERT(parameters.size() == 1, "Need exactly one parameter for proper dihedrals"); + mtop_->ffparams.iparams[tp].harmonic.krA = + mtop_->ffparams.iparams[tp].harmonic.krB = + ltop_->idef.iparams[tp].harmonic.krA = + ltop_->idef.iparams[tp].harmonic.krB = + gmx::doubleFromString(parameters[0].c_str()); + } + break; + default: + gmx_fatal(FARGS, "Unsupported dihedral type %s in UpdateIdef", + fs->function().c_str()); + } } } } - } - break; - case eitVDW: - case eitVSITE2: - case eitVSITE3FAD: - case eitVSITE3OUT: - case eitCONSTR: - case eitNR: - default: break; + case eitLJ14: + case eitPOLARIZATION: + case eitVDW: + case eitVSITE2: + case eitVSITE3FAD: + case eitVSITE3OUT: + case eitCONSTR: + case eitNR: + default: + break; + } } } diff --git a/src/programs/alexandria/mymol_low.cpp b/src/programs/alexandria/mymol_low.cpp index 14cf4cb7f0..7bdce29d10 100644 --- a/src/programs/alexandria/mymol_low.cpp +++ b/src/programs/alexandria/mymol_low.cpp @@ -727,6 +727,26 @@ void nonbondedFromPdToMtop(gmx_mtop_t *mtop, } } +void polarizabilityFromPdToMtop(gmx_mtop_t *mtop, + gmx_localtop_t *ltop, + t_atoms *atoms, + const Poldata *pd, + unsigned int ft) +{ + for (auto i = 0; i < ltop->idef.il[ft].nr; i += interaction_function[ft].nratoms+1) + { + auto tp = ltop->idef.il[ft].iatoms[i]; + auto ai = ltop->idef.il[ft].iatoms[i+1]; + auto polarUnit = string2unit(pd->getPolarUnit().c_str()); + double alpha, sigma; + if (pd->getAtypePol(*atoms->atomtype[ai], &alpha, &sigma)) + { + mtop->ffparams.iparams[tp].polarize.alpha = + ltop->idef.iparams[tp].polarize.alpha = convert2gmx(alpha, polarUnit); + } + } +} + void plist_to_mtop(const std::vector &plist, gmx_mtop_t *mtop_) { diff --git a/src/programs/alexandria/mymol_low.h b/src/programs/alexandria/mymol_low.h index cb07e870fa..4c52f8aeea 100644 --- a/src/programs/alexandria/mymol_low.h +++ b/src/programs/alexandria/mymol_low.h @@ -39,6 +39,7 @@ #include #include "gromacs/mdlib/mdatoms.h" +#include "gromacs/mdlib/shellfc.h" #include "gromacs/mdtypes/forcerec.h" #include "gromacs/pbcutil/pbc.h" @@ -126,6 +127,13 @@ void nonbondedFromPdToMtop(gmx_mtop_t *mtop, t_atoms *atoms, const Poldata *pd, t_forcerec *fr); + + +void polarizabilityFromPdToMtop(gmx_mtop_t *mtop, + gmx_localtop_t *ltop, + t_atoms *atoms, + const Poldata *pd, + unsigned int ft); void plist_to_mtop(const std::vector &plist, gmx_mtop_t *mtop_); diff --git a/src/programs/alexandria/poldata.cpp b/src/programs/alexandria/poldata.cpp index b698dfa457..6797598699 100644 --- a/src/programs/alexandria/poldata.cpp +++ b/src/programs/alexandria/poldata.cpp @@ -206,6 +206,26 @@ bool Poldata::atypeToPtype(const std::string &atype, return false; } +bool Poldata::ztypeToPtype(const std::string &ztype, + std::string &ptype) const +{ + if (ztype.size() == 0) + { + return false; + } + auto ai = std::find_if(alexandria_.begin(), alexandria_.end(), + [ztype](Ffatype const &fa) + { + return fa.getZtype().compare(ztype) == 0; + }); + if (ai != alexandria_.end() && ai->getPtype().size() > 0) + { + ptype = ai->getPtype(); + return true; + } + return false; +} + bool Poldata::getAtypePol(const std::string &atype, double *polar, double *sigPol) const @@ -218,6 +238,18 @@ bool Poldata::getAtypePol(const std::string &atype, return false; } +bool Poldata::getZtypePol(const std::string &ztype, + double *polar, + double *sigPol) const +{ + std::string ptype; + if (ztypeToPtype(ztype, ptype)) + { + return getPtypePol(ptype, polar, sigPol); + } + return false; +} + bool Poldata::getPtypePol(const std::string &ptype, double *polar, double *sigPol) const @@ -1080,6 +1112,71 @@ void Poldata::broadcast_eemprop(const t_commrec *cr) } } +void Poldata::broadcast_ptype(const t_commrec *cr) +{ + size_t nptype; + const int src = 0; + if (MASTER(cr)) + { + for (auto dest = 1; dest < cr->nnodes; dest++) + { + auto cs = gmx_send_data(cr, dest); + if (CS_OK == cs) + { + if (nullptr != debug) + { + fprintf(debug, "Going to update Poldata::ptype on node %d\n", dest); + } + /*Send Eemprops*/ + gmx_send_int(cr, dest, ptype_.size()); + for (auto &ptype : ptype_) + { + cs = ptype.Send(cr, dest); + } + } + gmx_send_done(cr, dest); + } + } + else + { + auto cs = gmx_recv_data(cr, src); + if (CS_OK == cs) + { + /*Receive Eemprops*/ + nptype = gmx_recv_int(cr, src); + ptype_.clear(); + for (size_t n = 0; (CS_OK == cs) && (n < nptype); n++) + { + Ptype ptype; + cs = ptype.Receive(cr, src); + if (CS_OK == cs) + { + ptype_.push_back(ptype); + if (nullptr != debug) + { + fprintf(debug, "Poldata::ptype is updated on node %d\n", cr->nodeid); + } + } + else + { + if (nullptr != debug) + { + fprintf(debug, "Could not update Poldata::ptype on node %d\n", cr->nodeid); + } + } + } + } + else + { + if (nullptr != debug) + { + fprintf(debug, "Could not update Poldata::ptype on node %d\n", cr->nodeid); + } + } + gmx_recv_data(cr, src); + } +} + void Poldata::broadcast(const t_commrec *cr) { const int src = 0; diff --git a/src/programs/alexandria/poldata.h b/src/programs/alexandria/poldata.h index b0ab1b1293..115f085c27 100644 --- a/src/programs/alexandria/poldata.h +++ b/src/programs/alexandria/poldata.h @@ -384,6 +384,14 @@ class Poldata bool atypeToPtype(const std::string &atype, std::string &ptype) const; + /*! \brief + * Return the poltype corresponding to ztype and true if successful + * + * \param[in] ztype Zeta type + * \param[out] ptype Polarizability type. + */ + bool ztypeToPtype(const std::string &ztype, + std::string &ptype) const; /*! \brief * Return the bond type corresponding to atom type and true if successful @@ -397,10 +405,16 @@ class Poldata /* Return 1 if OK, 0 if not found */ bool getPtypePol(const std::string &ptype, - double *polarizability, double *sigPol) const; + double *polarizability, + double *sigPol) const; bool getAtypePol(const std::string &atype, - double *polarizability, double *sigPol) const; + double *polarizability, + double *sigPol) const; + + bool getZtypePol(const std::string &ztype, + double *polarizability, + double *sigPol) const; void addMiller(const std::string &miller, int atomnumber, @@ -723,6 +737,9 @@ class Poldata //! Spread eemprop from master to slave nodes void broadcast_eemprop(const t_commrec *cr); + + //! Spread ptype from master to slave nodes + void broadcast_ptype(const t_commrec *cr); CommunicationStatus Send(const t_commrec *cr, int dest); diff --git a/src/programs/alexandria/tune_eem.cpp b/src/programs/alexandria/tune_eem.cpp index 90ba2e4f9e..1afe3967fb 100644 --- a/src/programs/alexandria/tune_eem.cpp +++ b/src/programs/alexandria/tune_eem.cpp @@ -310,6 +310,10 @@ double OptACM::calcDeviation() if (PAR(commrec()) && !final()) { poldata()->broadcast_eemprop(commrec()); + if (bFitAlpha_) + { + poldata()->broadcast_ptype(commrec()); + } } resetEnergies(); for (auto &mymol : mymols()) @@ -542,9 +546,9 @@ void OptACM::polData2TuneACM(real factor) } if (bFitAlpha_) { - auto alpha = 0.0; - auto sigma = 0.0; - if (poldata()->getAtypePol(ai->name(), &alpha, &sigma)) + auto alpha = 0.0; + auto sigma = 0.0; + if (poldata()->getZtypePol(ai->name(), &alpha, &sigma)) { if (0 != alpha) { @@ -556,6 +560,10 @@ void OptACM::polData2TuneACM(real factor) gmx_fatal(FARGS, "Polarizability is zero for atom %s\n", ai->name().c_str()); } } + else + { + gmx_fatal(FARGS, "No Ptype for zeta type %s\n", ai->name().c_str()); + } } } } @@ -624,14 +632,14 @@ void OptACM::toPolData(const std::vector &changed) if (bFitAlpha_) { std::string ptype; - if (pd->atypeToPtype(ai->name(), ptype)) + if (pd->ztypeToPtype(ai->name(), ptype)) { pd->setPtypePolarizability(ptype, param[n], psigma[n]); n++; } else { - gmx_fatal(FARGS, "No Ptype for atom type %s\n", ai->name().c_str()); + gmx_fatal(FARGS, "No Ptype for zeta type %s\n", ai->name().c_str()); } } } diff --git a/src/programs/alexandria/tune_zeta.cpp b/src/programs/alexandria/tune_zeta.cpp index c90a4f8d62..f4069e5e58 100644 --- a/src/programs/alexandria/tune_zeta.cpp +++ b/src/programs/alexandria/tune_zeta.cpp @@ -296,7 +296,11 @@ double OptZeta::calcDeviation() } if (PAR(commrec()) && !final()) { - poldata()->broadcast(commrec()); + poldata()->broadcast_eemprop(commrec()); + if (bFitAlpha_) + { + poldata()->broadcast_ptype(commrec()); + } } resetEnergies(); for (auto &mymol : mymols()) -- 2.11.4.GIT