Inniag / openmm-scripts-amoeba

:microscope: OpenMM scripts for polarisable molecular dynamics with the AMOEBA force field
MIT License
23 stars 4 forks source link

AMOEBA for small molecule from Tinker poltype2 #4

Open pradipchm opened 1 month ago

pradipchm commented 1 month ago

Hi, First of all thanks for providing the thinker2openmm code for amoeba force field. I am very new to the amoeba force field. I am using your code to covert the poltype2 format to openmm xml using your code. I changed few lines to align to the poltype2 format. I have doubt about few things. I would appreciate your suggestions and help. In the mutlipole section, the number of parameters can vary, like 13 or 14 or 15, based on that I have changed the code.

# [406] = [[7]]
multipole   406  405                    0.50670
                                        0.00000    0.00000    0.13788
                                        0.55354
                                        0.00000    0.55354
                                        0.00000    0.00000   -1.10708
# [407] = [[8]]
multipole   407  417  401              -0.14418
                                        0.03235    0.00000    0.15851
                                       -0.12079
                                        0.00000   -0.16191
                                        0.05107    0.00000    0.28270

the new function looks like

def parse_prm_multipole(self, split_line, f):

        # Get the length of split_line and adjust variables accordingly
        atom_type = split_line[1]

        if len(split_line) == 4:
            # Case: "multipole 410 409 "
            kz = split_line[2]
            c_0 = split_line[3]
            kx = None  # No kx for this case
            ky = None  # No ky for this case

        elif len(split_line) == 5:
            # Case: "multipole 410 409 411"
            kz = split_line[2]
            kx = split_line[3]
            c_0 = split_line[4]
            ky = None  # No ky for this case

        elif len(split_line) == 6:
            # Case: "multipole 410 409 411 412"
            kz = split_line[2]
            kx = split_line[3]
            ky = split_line[4]
            c_0 = split_line[5]

        else:
            raise ValueError(f"Unexpected number of parameters for multipole: {split_line}")
        # Get second line parameters (dipole moment):
        line = next(f).strip()  # read next line for dipole
        split_line = shlex.split(line)
        d_1 = split_line[0]
        d_2 = split_line[1]
        d_3 = split_line[2]

        # Get third line parameters (quadrupole moments):
        line = next(f).strip()  # read next line for quadrupole
        split_line = shlex.split(line)
        q_11 = split_line[0]

        # Get fourth line parameters (quadrupole moments):
        line = next(f).strip()  # read next line for quadrupole
        split_line = shlex.split(line)
        q_21 = split_line[0]
        q_22 = split_line[1]

        # Get fifth line parameters (quadrupole moments):
        line = next(f).strip()  # read next line for quadrupole
        split_line = shlex.split(line)
        q_31 = split_line[0]
        q_32 = split_line[1]
        q_33 = split_line[2]
        # create multipole record:
          self.prm_multipole_records.append(TinkerPrmMultipole(
            atom_type=atom_type,
            kx=kx,
            kz=kz,
            ky=ky,
            c_0=c_0,
            d_1=d_1,
            d_2=d_2,
            d_3=d_3,
            q_11=q_11,
            q_21=q_21,
            q_22=q_22,
            q_31=q_31,
            q_32=q_32,
            q_33=q_33))  

Then while reading the prm file there are conditions that I don't have in my file, what are these parameters signifies, like

# parse global van der Waals record types:
                elif split_line[0] == "vdwtype":
                    self.vdwtype = split_line[1]
                elif split_line[0] == "radiusrule":
                    self.radiusrule = split_line[1]
                elif split_line[0] == "radiustype":
                    self.radiustype = split_line[1]
                elif split_line[0] == "radiussize":
                    self.radiussize = split_line[1]
                elif split_line[0] == "epsilonrule":
                    self.epsilonrule = split_line[1]
                elif split_line[0] == "vdw-12-scale":
                    self.vdw_12_scale = split_line[1]
                elif split_line[0] == "vdw-13-scale":
                    self.vdw_13_scale = split_line[1]
                elif split_line[0] == "vdw-14-scale":
                    self.vdw_14_scale = split_line[1]
                elif split_line[0] == "vdw-15-scale":
                    self.vdw_15_scale = split_line[1]
                # parse global bond force records:
                elif split_line[0] == "bond-cubic":
                    self.bond_cubic = split_line[1]
                elif split_line[0] == "bond-quartic":
                    self.bond_quartic = split_line[1]

                # parse global angle force records:
                elif split_line[0] == "angle-cubic":
                    self.angle_cubic = split_line[1]
                elif split_line[0] == "angle-quartic":
                    self.angle_quartic = split_line[1]
                elif split_line[0] == "angle-pentic":
                    self.angle_pentic = split_line[1]
                elif split_line[0] == "angle-sextic":
                    self.angle_sextic = split_line[1]

                # parse out-of-plane-bending force parameters:
                elif split_line[0] == "opbendtype":
                    self.opbend_type = split_line[1]
                elif split_line[0] == "opbend-cubic":
                    self.opbend_cubic = split_line[1]
                elif split_line[0] == "opbend-quartic":
                    self.opbend_quartic = split_line[1]
                elif split_line[0] == "opbend-pentic":
                    self.opbend_pentic = split_line[1]
                elif split_line[0] == "opbend-sextic":
                    self.opbend_sextic = split_line[1]

                # parse global torsion force records:
                elif split_line[0] == "torsionunit":
                    self.torsion_unit = split_line[1]

                # parse global mutlipole parameters:
                elif split_line[0] == "direct-11-scale":
                    self.direct_11_scale = split_line[1]
                elif split_line[0] == "direct-12-scale":
                    self.direct_12_scale = split_line[1]
                elif split_line[0] == "direct-13-scale":
                    self.direct_13_scale = split_line[1]
                elif split_line[0] == "direct-14-scale":
                    self.direct_14_scale = split_line[1]
                elif split_line[0] == "mpole-12-scale":
                    self.mpole_12_scale = split_line[1]
                elif split_line[0] == "mpole-13-scale":
                    self.mpole_13_scale = split_line[1]
                elif split_line[0] == "mpole-14-scale":
                    self.mpole_14_scale = split_line[1]
                elif split_line[0] == "mpole-15-scale":
                    self.mpole_15_scale = split_line[1]
                elif split_line[0] == "mutual-11-scale":
                    self.mutual_11_scale = split_line[1]
                elif split_line[0] == "mutual-12-scale":
                    self.mutual_12_scale = split_line[1]
                elif split_line[0] == "mutual-13-scale":
                    self.mutual_13_scale = split_line[1]
                elif split_line[0] == "mutual-14-scale":
                    self.mutual_14_scale = split_line[1]
               elif split_line[0] == "polar-15-scale":
                    self.polar_15_scale = split_line[1]
                elif split_line[0] == "polar-14-intra":
                    self.polar_14_intra = split_line[1]

Are these parameters important. If you would like to share me any documentation or information about openmm amoeba that might help me to understand.