diff --git a/corelib/src/libs/SireIO/grotop.cpp b/corelib/src/libs/SireIO/grotop.cpp index 43c169a15..87908ea6b 100644 --- a/corelib/src/libs/SireIO/grotop.cpp +++ b/corelib/src/libs/SireIO/grotop.cpp @@ -2869,6 +2869,7 @@ static QStringList writeAtomTypes(QMap, GroMolType> &moltyps QString particle_type = "A"; // A is for Atom + // This is a dummy atom. if (elem.nProtons() == 0 and lj.isDummy()) { if (is_perturbable) @@ -2877,6 +2878,9 @@ static QStringList writeAtomTypes(QMap, GroMolType> &moltyps // Only label dummies for regular simulations. else if (not was_perturbable) particle_type = "D"; + + // Flag that we need to update the atoms. + update_atoms0 = true; } // This is a new atom type. @@ -2893,6 +2897,15 @@ static QStringList writeAtomTypes(QMap, GroMolType> &moltyps // Hash the atom type against its parameter string, minus the type. param_hash.insert(atomtypes[atomtype].mid(6), atomtype); + + if (update_atoms0) + { + // Set the type. + atom.setAtomType(atomtype); + + // Update the atoms in the vector. + atoms[i] = atom; + } } // This type has been seen before. else @@ -2977,6 +2990,17 @@ static QStringList writeAtomTypes(QMap, GroMolType> &moltyps update_atoms0 = true; } } + else + { + if (update_atoms0) + { + // Set the type. + atom.setAtomType(atomtype); + + // Update the atoms in the vector. + atoms[i] = atom; + } + } } } @@ -3019,9 +3043,15 @@ static QStringList writeAtomTypes(QMap, GroMolType> &moltyps QString particle_type = "A"; // A is for Atom + // This is a dummy atom. if (elem.nProtons() == 0 and lj.isDummy()) + { atomtype += "_du"; + // Flag that we need to update the atoms. + update_atoms1 = true; + } + // This is a new atom type. if (not atomtypes.contains(atomtype)) { @@ -3036,6 +3066,15 @@ static QStringList writeAtomTypes(QMap, GroMolType> &moltyps // Hash the atom type against its parameter string, minus the type. param_hash.insert(atomtypes[atomtype].mid(6), atomtype); + + if (update_atoms1) + { + // Set the type. + atom.setAtomType(atomtype); + + // Update the atoms in the vector. + atoms[i] = atom; + } } // This type has been seen before. @@ -3121,6 +3160,17 @@ static QStringList writeAtomTypes(QMap, GroMolType> &moltyps update_atoms1 = true; } } + else + { + if (update_atoms1) + { + // Set the type. + atom.setAtomType(atomtype); + + // Update the atoms in the vector. + atoms[i] = atom; + } + } } } @@ -3217,14 +3267,6 @@ static QStringList writeMolType(const QString &name, const GroMolType &moltype, elem1 = Element::elementWithMass(mol.property("mass1").asA()[cgatomidx]); } - // Update the atom types. - - if (elem0.nProtons() == 0) - atomtype0 += "_du"; - - if (elem1.nProtons() == 0) - atomtype1 += "_du"; - QString resnum = QString::number(atom0.residueNumber().value()); if (not atom0.chainName().isNull()) diff --git a/corelib/src/libs/SireMove/openmmpmefep.cpp b/corelib/src/libs/SireMove/openmmpmefep.cpp index 6c2c8096d..1332bf8ae 100644 --- a/corelib/src/libs/SireMove/openmmpmefep.cpp +++ b/corelib/src/libs/SireMove/openmmpmefep.cpp @@ -836,15 +836,80 @@ void OpenMMPMEFEP::initialise(bool fullPME) qDebug() << "\n\nRestraint is ON\n\n"; } + /************************************RECEPTOR-lIGAND RESTRAINTS**************************/ + // Check if we are in turn on receptor-ligand restraint mode + + bool turn_on_restraints_mode{false}; + + for (int i = 0; i < nmols; i++) + { + Molecule molecule = moleculegroup.moleculeAt(i).molecule(); + + if (molecule.hasProperty("turn_on_restraints_mode")) + { + turn_on_restraints_mode = true; // Lambda will be used to turn on the receptor-ligand restraints + if (Debug) + qDebug() << "Lambda will be used to turn on the receptor-ligand restraints"; + break; // We've found the solute - exit loop over molecules in system. + } + } + /*** BOND LINK FORCE FIELD ***/ - /* NOTE: CustomBondForce does not (OpenMM 6.2) apply PBC checks so code will be buggy if - restraints involve one atom that diffuses out of the box. */ + /* FC 12/21 CustomBondForce now (OpenMM 7.4.0) allows application of PBC checks*/ - auto custom_link_bond = new OpenMM::CustomBondForce("kl * max(0, d - dl*dl);" - "d = (r-reql) * (r-reql)"); + OpenMM::CustomBondForce * custom_link_bond = new OpenMM::CustomBondForce("delta(min(0, r_eff))*(lamrest^5)*kl*r_eff^2;" + "r_eff=abs(r-reql)-dl"); custom_link_bond->addPerBondParameter("reql"); custom_link_bond->addPerBondParameter("kl"); custom_link_bond->addPerBondParameter("dl"); + custom_link_bond->setUsesPeriodicBoundaryConditions(true); + // If in turn on receptor-ligand restraints mode, default value of lamrest needs to be lambda, because + // the default value is used for the first nrg_freq timesteps before being set by updateOpenMMContextLambda + if (turn_on_restraints_mode) + custom_link_bond->addGlobalParameter("lamrest", current_lambda); + // We are not in turn on receptor-ligand restraints mode - set lamrest to 1 + else + custom_link_bond->addGlobalParameter("lamrest", 1); + + /****************************************BORESCH DISTANCE POTENTIAL*****************************/ + + OpenMM::CustomBondForce *custom_boresch_dist_rest = + new OpenMM::CustomBondForce("lamrest*force_const*(r-equil_val)^2"); + custom_boresch_dist_rest->addPerBondParameter("force_const"); + custom_boresch_dist_rest->addPerBondParameter("equil_val"); + custom_boresch_dist_rest->setUsesPeriodicBoundaryConditions(true); + if (turn_on_restraints_mode) + custom_boresch_dist_rest->addGlobalParameter("lamrest", current_lambda); + // We are not in turn on receptor-ligand restraints mode - set lamrest to 1 + else + custom_boresch_dist_rest->addGlobalParameter("lamrest", 1); + + /****************************************BORESCH ANGLE POTENTIAL*****************************/ + + OpenMM::CustomAngleForce *custom_boresch_angle_rest = + new OpenMM::CustomAngleForce("lamrest*force_const*(theta-equil_val)^2"); + custom_boresch_angle_rest->addPerAngleParameter("force_const"); + custom_boresch_angle_rest->addPerAngleParameter("equil_val"); + custom_boresch_angle_rest->setUsesPeriodicBoundaryConditions(true); + if (turn_on_restraints_mode) + custom_boresch_angle_rest->addGlobalParameter("lamrest", current_lambda); + // We are not in turn on receptor-ligand restraints mode - set lamrest to 1 + else + custom_boresch_angle_rest->addGlobalParameter("lamrest", 1); + + /****************************************BORESCH DIHEDRAL POTENTIAL*****************************/ + + OpenMM::CustomTorsionForce *custom_boresch_dihedral_rest = + new OpenMM::CustomTorsionForce("lamrest*force_const*min(dtheta, 2*pi-dtheta)^2;" + "dtheta = abs(theta-equil_val); pi = 3.1415926535"); + custom_boresch_dihedral_rest->addPerTorsionParameter("force_const"); + custom_boresch_dihedral_rest->addPerTorsionParameter("equil_val"); + custom_boresch_dihedral_rest->setUsesPeriodicBoundaryConditions(true); + if (turn_on_restraints_mode) + custom_boresch_dihedral_rest->addGlobalParameter("lamrest", current_lambda); + // We are not in turn on receptor-ligand restraints mode - set lamrest to 1 + else + custom_boresch_dihedral_rest->addGlobalParameter("lamrest", 1); /*** BUILD OpenMM SYSTEM ***/ @@ -970,7 +1035,10 @@ void OpenMMPMEFEP::initialise(bool fullPME) // Molecule solutemol = solute.moleculeAt(0).molecule(); int nions = 0; - QVector perturbed_energies_tmp{false, false, false, false, false, false, false, false, false}; + QVector perturbed_energies_tmp(10); + + for (int i = 0; i < perturbed_energies_tmp.size(); i++) + perturbed_energies_tmp[i] = false; // the default AMBER 1-4 scaling factors double const Coulomb14Scale = 1.0 / 1.2; @@ -1853,7 +1921,7 @@ void OpenMMPMEFEP::initialise(bool fullPME) dihedral_pert_list.append(DihedralID(four.atom0(), four.atom1(), four.atom2(), four.atom3())); dihedral_pert_swap_list.append( - DihedralID(four.atom3(), four.atom1(), four.atom2(), four.atom0())); + DihedralID(four.atom3(), four.atom2(), four.atom1(), four.atom0())); improper_pert_list.append(ImproperID(four.atom0(), four.atom1(), four.atom2(), four.atom3())); improper_pert_swap_list.append( @@ -2501,6 +2569,14 @@ void OpenMMPMEFEP::initialise(bool fullPME) } } // if (!fullPME) + + if (turn_on_restraints_mode) + { + perturbed_energies_tmp[9] = true; //Lambda will be used to turn on the receptor-ligand restraints + if (Debug) + qDebug() << "Added Perturbed Receptor-Ligand Restraint energy term"; + } + perturbed_energies = perturbed_energies_tmp; // IMPORTANT: PERTURBED ENERGY TORSIONS ARE ADDED ABOVE @@ -2556,6 +2632,174 @@ void OpenMMPMEFEP::initialise(bool fullPME) } // end of bond link flag + bool UseBoresch_flag = true; + + // Boresch Restraints. All the information is stored in the solute only. + + if (UseBoresch_flag == true) + { + bool found_solute{false}; + for (int i = 0; i < nmols; i++) + { + Molecule molecule = moleculegroup.moleculeAt(i).molecule(); + + bool has_boresch_dist = molecule.hasProperty("boresch_dist_restraint"); + bool has_boresch_angle = molecule.hasProperty("boresch_angle_restraints"); + bool has_boresch_dihedral = molecule.hasProperty("boresch_dihedral_restraints"); + + if (has_boresch_dist) + { + found_solute = true; // We have found the solute, but before breaking we must also check + // if there are Boresch angle and torsion restraints. + + if (Debug) + { + qDebug() << "Boresch distance restraint properties stored = true"; + qDebug() << "Boresch angle restraint properties stored = " << has_boresch_angle; + qDebug() << "Boresch dihedral restraint properties stored = " << has_boresch_dihedral; + } + + std::vector custom_boresch_dist_par(2); + + const auto boresch_dist_prop = molecule.property("boresch_dist_restraint").asA(); + + const auto atomnum0 = boresch_dist_prop.property(QString("AtomNum0")).asA().toInt(); + const auto atomnum1 = boresch_dist_prop.property(QString("AtomNum1")).asA().toInt(); + const auto force_const = + boresch_dist_prop.property(QString("force_const")).asA().toDouble(); + const auto equil_val = + boresch_dist_prop.property(QString("equil_val")).asA().toDouble(); + + const auto openmmindex0 = AtomNumToOpenMMIndex[atomnum0]; + const auto openmmindex1 = AtomNumToOpenMMIndex[atomnum1]; + + custom_boresch_dist_par[0] = + force_const * (OpenMM::KJPerKcal * OpenMM::AngstromsPerNm * OpenMM::AngstromsPerNm); // force_const + custom_boresch_dist_par[1] = equil_val * OpenMM::NmPerAngstrom; // equil_val + + if (Debug) + { + qDebug() << "Boresch distance restraint implemented"; + qDebug() << "atomnum0 = " << atomnum0 << " openmmindex0 =" << openmmindex0; + qDebug() << "atomnum1 = " << atomnum1 << " openmmindex1 =" << openmmindex1; + qDebug() << "force_const = " << force_const << " equil_val = " << equil_val; + } + + custom_boresch_dist_rest->addBond(openmmindex0, openmmindex1, custom_boresch_dist_par); + + system_openmm->addForce(custom_boresch_dist_rest); + } + + if (has_boresch_angle) + { + std::vector custom_boresch_angle_par(2); + + const auto boresch_angle_prop = molecule.property("boresch_angle_restraints").asA(); + + const auto n_angles = + boresch_angle_prop.property(QString("n_boresch_angle_restraints")).asA().toInt(); + + if (Debug) + qDebug() << "Number of Boresch angle restraints = " << n_angles; + + for (int i = 0; i < n_angles; i++) + { + const auto atomnum0 = + boresch_angle_prop.property(QString("AtomNum0-%1").arg(i)).asA().toInt(); + const auto atomnum1 = + boresch_angle_prop.property(QString("AtomNum1-%1").arg(i)).asA().toInt(); + const auto atomnum2 = + boresch_angle_prop.property(QString("AtomNum2-%1").arg(i)).asA().toInt(); + const auto force_const = + boresch_angle_prop.property(QString("force_const-%1").arg(i)).asA().toDouble(); + const auto equil_val = + boresch_angle_prop.property(QString("equil_val-%1").arg(i)).asA().toDouble(); + + const auto openmmindex0 = AtomNumToOpenMMIndex[atomnum0]; + const auto openmmindex1 = AtomNumToOpenMMIndex[atomnum1]; + const auto openmmindex2 = AtomNumToOpenMMIndex[atomnum2]; + + custom_boresch_angle_par[0] = force_const * (OpenMM::KJPerKcal); // force_const + custom_boresch_angle_par[1] = equil_val; // equil_val + + if (Debug) + { + qDebug() << "atomnum0 = " << atomnum0 << " openmmindex0 =" << openmmindex0; + qDebug() << "atomnum1 = " << atomnum1 << " openmmindex1 =" << openmmindex1; + qDebug() << "atomnum2 = " << atomnum2 << " openmmindex2 =" << openmmindex2; + qDebug() << "force_const = " << force_const << " equil_val = " << equil_val; + } + + custom_boresch_angle_rest->addAngle(openmmindex0, openmmindex1, openmmindex2, + custom_boresch_angle_par); + } + + system_openmm->addForce(custom_boresch_angle_rest); + } + + if (has_boresch_dihedral) + { + std::vector custom_boresch_dihedral_par(2); + + const auto boresch_dihedral_prop = molecule.property("boresch_dihedral_restraints").asA(); + + const auto n_dihedrals = boresch_dihedral_prop.property(QString("n_boresch_dihedral_restraints")) + .asA() + .toInt(); + + if (Debug) + qDebug() << "Number of Boresch dihedral restraints = " << n_dihedrals; + + for (int i = 0; i < n_dihedrals; i++) + { + const auto atomnum0 = + boresch_dihedral_prop.property(QString("AtomNum0-%1").arg(i)).asA().toInt(); + const auto atomnum1 = + boresch_dihedral_prop.property(QString("AtomNum1-%1").arg(i)).asA().toInt(); + const auto atomnum2 = + boresch_dihedral_prop.property(QString("AtomNum2-%1").arg(i)).asA().toInt(); + const auto atomnum3 = + boresch_dihedral_prop.property(QString("AtomNum3-%1").arg(i)).asA().toInt(); + const auto force_const = boresch_dihedral_prop.property(QString("force_const-%1").arg(i)) + .asA() + .toDouble(); + const auto equil_val = boresch_dihedral_prop.property(QString("equil_val-%1").arg(i)) + .asA() + .toDouble(); + + const auto openmmindex0 = AtomNumToOpenMMIndex[atomnum0]; + const auto openmmindex1 = AtomNumToOpenMMIndex[atomnum1]; + const auto openmmindex2 = AtomNumToOpenMMIndex[atomnum2]; + const auto openmmindex3 = AtomNumToOpenMMIndex[atomnum3]; + + custom_boresch_dihedral_par[0] = force_const * (OpenMM::KJPerKcal); // force_const + custom_boresch_dihedral_par[1] = equil_val; // equil_val + + if (Debug) + { + qDebug() << "atomnum0 = " << atomnum0 << " openmmindex0 =" << openmmindex0; + qDebug() << "atomnum1 = " << atomnum1 << " openmmindex1 =" << openmmindex1; + qDebug() << "atomnum2 = " << atomnum2 << " openmmindex2 =" << openmmindex2; + qDebug() << "atomnum3 = " << atomnum3 << " openmmindex3 =" << openmmindex3; + qDebug() << "force_const = " << force_const << " equil_val = " << equil_val; + } + + custom_boresch_dihedral_rest->addTorsion(openmmindex0, openmmindex1, openmmindex2, openmmindex3, + custom_boresch_dihedral_par); + } + + system_openmm->addForce(custom_boresch_dihedral_rest); + } + + if (found_solute) + break; // We've found the molecule, exit the outer loop. If a molecule has Boresch + // distance restraints it must be the solute, but we cannot break immediately + // because it may also have angle/ dihedral restraints + + } // End of loop over molecules in system + + } // End of Boresch flag + this->openmm_system = system_openmm; this->isSystemInitialised = true; } // OpenMMPMEFEP::initialise END @@ -3414,6 +3658,10 @@ void OpenMMPMEFEP::updateOpenMMContextLambda(double lambda) if (perturbed_energies[7]) openmm_context->setParameter("lamdih", lambda); // Torsions + // RECEPTOR-LIGAND RESTRAINTS + if (perturbed_energies[9]) + openmm_context->setParameter("lamrest", lambda); //Receptor-ligand restraints + // lambda for the offsets (linear scaling) of the charges in // reciprocal space openmm_context->setParameter("lambda_offset", lambda); diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index 3a1d98aba..915311990 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -23,6 +23,9 @@ organisation on `GitHub `__. * Set ``IFBOX`` pointer to 3 for general triclinic boxes in ``sire.IO.AmberPrm`` parser. * Only excluded nonbonded interactions between from_ghost and to_ghost atoms if they are in the same molecule. * Add Docker support for building wrappers on Linux x86. +* Add support for boresch restraints to PME. +* Port SOMD torsion fix to PME code. +* Fix issues with ``atomtype`` and ``atom`` records for dummy atoms in GROMACS topology files. `2024.2.0 `__ - June 2024 diff --git a/tests/io/test_grotop.py b/tests/io/test_grotop.py index 5a23a87e4..09d93a603 100644 --- a/tests/io/test_grotop.py +++ b/tests/io/test_grotop.py @@ -39,3 +39,131 @@ def test_posre(): # Make sure we can parse a file with BioSimSpace position restraint include # directives. mols = sr.load_test_files("posre.top") + + +def test_fep_atoms(): + """ + Test that GROMACS FEP atomtypes and atoms are created correctly. + """ + + import os + import tempfile + + atomtypes = """[ atomtypes ] + ; name at.num mass charge ptype sigma epsilon + C1 6 12.010700 0.000000 A 0.348065 0.363503 + C1_du 0 0.000000 0.000000 A 0.348065 0.000000 + C2 6 12.010700 0.000000 A 0.337953 0.455389 + H1 1 1.007940 0.000000 A 0.110343 0.058956 + H1_du 0 0.000000 0.000000 A 0.110343 0.000000 + H2 1 1.007940 0.000000 A 0.257258 0.065318 + H2_du 0 0.000000 0.000000 A 0.264454 0.000000 + H2_dux 0 0.000000 0.000000 A 0.257258 0.000000 + H3 1 1.007940 0.000000 A 0.264454 0.066021 + H3_du 0 0.000000 0.000000 A 0.258323 0.000000 + H4 1 1.007940 0.000000 A 0.258323 0.068656 + H5 1 1.007940 0.000000 A 0.245363 0.054840 + N1 7 14.006700 0.000000 A 0.320688 0.701621 + O1 8 15.999400 0.000000 A 0.303981 0.879502 + O1_du 0 0.000000 0.000000 A 0.303981 0.000000 + O2 8 15.999400 0.000000 A 0.302511 0.704858 + """ + + atoms = """[ atoms ] + ; nr type0 resnr residue atom cgnr charge0 mass0 type1 charge1 mass1 + 1 O1 1 LIG O1x 1 -0.803190 15.999430 O1_du 0.000000 15.999430 + 2 C1 1 LIG C1x 2 0.916780 12.010780 N1 -0.662150 14.006720 + 3 O1 1 LIG O2x 3 -0.803190 15.999430 O1_du 0.000000 15.999430 + 4 C1 1 LIG C2x 4 0.082410 12.010780 C1 -0.105710 12.010780 + 5 N1 1 LIG N1x 5 -0.564860 14.006720 N1 -0.484790 14.006720 + 6 H1 1 LIG H1x 6 0.449360 1.007947 H1 0.407460 1.007947 + 7 C1 1 LIG C3x 7 0.061040 12.010780 C1 0.001920 12.010780 + 8 C1 1 LIG C4x 8 -0.087870 12.010780 C1 -0.087450 12.010780 + 9 C1 1 LIG C5x 9 -0.101450 12.010780 N1 -0.015090 14.006720 + 10 C1 1 LIG C6x 10 -0.165020 12.010780 C1 -0.045730 12.010780 + 11 H2 1 LIG H2x 11 0.106400 1.007947 H5 0.196160 1.007947 + 12 C1 1 LIG C7x 12 -0.107040 12.010780 C1_du 0.000000 12.010780 + 13 H2 1 LIG H3x 13 0.120980 1.007947 H2_dux 0.000000 1.007947 + 14 C1 1 LIG C8x 14 -0.057540 12.010780 C1 -0.204760 12.010780 + 15 C1 1 LIG C9x 15 -0.203310 12.010780 C1 0.114810 12.010780 + 16 C2 1 LIG C10x 16 0.013180 12.010780 C2 -0.054450 12.010780 + 17 H3 1 LIG H4x 17 0.044890 1.007947 H3 0.058540 1.007947 + 18 H3 1 LIG H5x 18 0.044890 1.007947 H3 0.058540 1.007947 + 19 C2 1 LIG C11x 19 -0.084960 12.010780 C2 -0.080770 12.010780 + 20 H3 1 LIG H6x 20 0.058790 1.007947 H3 0.071370 1.007947 + 21 H3 1 LIG H7x 21 0.058790 1.007947 H3 0.071370 1.007947 + 22 C2 1 LIG C12x 22 0.126550 12.010780 C2 0.130690 12.010780 + 23 H4 1 LIG H8x 23 0.037860 1.007947 H4 0.038330 1.007947 + 24 H4 1 LIG H9x 24 0.037860 1.007947 H4 0.038330 1.007947 + 25 O2 1 LIG O3x 25 -0.332540 15.999430 O2 -0.334080 15.999430 + 26 C1 1 LIG C13x 26 0.142880 12.010780 C1 0.108960 12.010780 + 27 C1 1 LIG C14x 27 -0.181170 12.010780 C1 -0.171770 12.010780 + 28 H2 1 LIG H10x 28 0.144890 1.007947 H2 0.138510 1.007947 + 29 C1 1 LIG C15x 29 -0.052890 12.010780 C1 -0.042000 12.010780 + 30 C2 1 LIG C16x 30 -0.051460 12.010780 C2 -0.059900 12.010780 + 31 H3 1 LIG H11x 31 0.040320 1.007947 H3 0.050460 1.007947 + 32 H3 1 LIG H12x 32 0.040320 1.007947 H3 0.050460 1.007947 + 33 H3 1 LIG H13x 33 0.040320 1.007947 H3 0.050460 1.007947 + 34 C1 1 LIG C17x 34 -0.175410 12.010780 C1 -0.147890 12.010780 + 35 H2 1 LIG H14x 35 0.122060 1.007947 H2 0.142460 1.007947 + 36 C1 1 LIG C18x 36 -0.101660 12.010780 C1 -0.094400 12.010780 + 37 H2 1 LIG H15x 37 0.123170 1.007947 H2 0.139870 1.007947 + 38 C1 1 LIG C19x 38 -0.184410 12.010780 C1 -0.174770 12.010780 + 39 H2 1 LIG H16x 39 0.144680 1.007947 H2 0.138300 1.007947 + 40 C2 1 LIG C20x 40 -0.038120 12.010780 C2 -0.032020 12.010780 + 41 H3 1 LIG H17x 41 0.030840 1.007947 H2_du 0.000000 1.007947 + 42 H3 1 LIG H18x 42 0.030840 1.007947 H2_du 0.000000 1.007947 + 43 H3 1 LIG H19x 43 0.030840 1.007947 H2_du 0.000000 1.007947 + 44 C2 1 LIG C21x 44 -0.035100 12.010780 C2 0.001470 12.010780 + 45 H3 1 LIG H20x 45 0.026750 1.007947 H2_du 0.000000 1.007947 + 46 H3 1 LIG H21x 46 0.026750 1.007947 H2_du 0.000000 1.007947 + 47 H3 1 LIG H22x 47 0.026750 1.007947 H2_du 0.000000 1.007947 + 48 H1_du 1 LIG H1x 48 0.000000 1.007947 H1 0.459340 1.007947 + 49 H1_du 1 LIG H2x 49 0.000000 1.007947 H1 0.459340 1.007947 + 50 H1_du 1 LIG H3x 50 0.000000 1.007947 H1 0.459340 1.007947 + 51 H2_du 1 LIG H19x 51 0.000000 1.007947 H3 0.062430 1.007947 + 52 H2_du 1 LIG H20x 52 0.000000 1.007947 H3 0.062430 1.007947 + 53 H2_du 1 LIG H21x 53 0.000000 1.007947 H3 0.062430 1.007947 + 54 H3_du 1 LIG H22x 54 0.000000 1.007947 H4 0.074650 1.007947 + 55 H3_du 1 LIG H23x 55 0.000000 1.007947 H4 0.074650 1.007947 + 56 H3_du 1 LIG H24x 56 0.000000 1.007947 H4 0.074650 1.007947 + """ + + # Load the molecule. + merged_mol = sr.load_test_files("merged_molecule_grotop.s3") + + # Save to a temporary file. + with tempfile.TemporaryDirectory() as tmpdir: + sr.save(merged_mol, os.path.join(tmpdir, "merged_molecule"), format="GroTop") + + # Read the [ atomtypes ] and [ atoms ] sections from the file. + with open(os.path.join(tmpdir, "merged_molecule.grotop"), "r") as f: + atomtypes_lines = [] + atoms_lines = [] + found_atomtypes = False + found_atoms = False + for line in f: + if line.startswith("[ atomtypes ]"): + found_atomtypes = True + elif line.startswith("[ atoms ]"): + found_atoms = True + + if found_atomtypes and not found_atoms: + if line != "\n": + atomtypes_lines.append(line) + else: + found_atomtypes = False + + if found_atoms: + if line != "\n": + atoms_lines.append(line) + else: + found_atoms = False + + # Check that the atomtypes and atoms are as expected. + atomtypes_list = atomtypes.split("\n") + for a, b in zip(atomtypes_list, atomtypes_lines): + assert a.strip() == b.strip() + atoms_list = atoms.split("\n") + for a, b in zip(atoms_list, atoms_lines): + assert a.strip() == b.strip()