openforcefield / openff-interchange

A project (and object) for storing, manipulating, and converting molecular mechanics data.
https://docs.openforcefield.org/projects/interchange
MIT License
69 stars 22 forks source link

Time value in inpcrd is messing up VMD visualization #854

Open GregorySchwing opened 10 months ago

GregorySchwing commented 10 months ago

Description to_inpcrd writes a time of 0.0

Reproduction

`
from openff.toolkit import ForceField, Molecule, Topology
from openff.units import Quantity, unit
from openff.interchange import Interchange
from openff.toolkit.typing.engines.smirnoff import ForceField
from openff.toolkit.typing.engines import smirnoff

def build_water(OH_bond_length=1, HH_bond_length=1.63298086):
    from math import cos, sin, acos, degrees
    calc_theta = lambda a, b: degrees(acos(1 - ((a * a) / (2 * b * b)))) if a != 0 else None    
    water_reference = Molecule.from_mapped_smiles("[O:1]([H:2])[H:3]")
    water_reference.atoms[0].name = "O"
    water_reference.atoms[1].name = "H1"
    water_reference.atoms[2].name = "H2"

    # Add ideal TIP5P geometry
    bond_length = Quantity(OH_bond_length, unit.angstrom)
    print("theta",calc_theta(HH_bond_length,OH_bond_length))
    print("rOH",bond_length)

    theta = Quantity(calc_theta(HH_bond_length,OH_bond_length), unit.degree).to(unit.radian)
    water_reference.add_conformer(
        bond_length
        * Quantity(
            [
                [0.0, 0.0, 0.0],
                [-sin(theta / 2), cos(theta / 2), 0.0],
                [sin(theta / 2), cos(theta / 2), 0.0],
            ]
        )
    )

    return water_reference

#print(smirnoff.get_available_force_fields())
print("tip3p-1.0.1.offxml")
water_interchange = ForceField("tip3p-1.0.1.offxml").create_interchange(Molecule.from_mapped_smiles(
    "[O:1]([H:2])[H:3]"
).to_topology())
OH_bond, HH_bond = water_interchange.collections['Constraints'].potentials.values()

# Build correct water molecule
water_reference = build_water(OH_bond.parameters['distance'],HH_bond.parameters['distance'])

# Scaffold interchange, with bond information needed to construct water
interchange = ForceField("openff-2.0.0.offxml").create_interchange(water_reference.to_topology())

# Assigns desired forcefield parameters to interchange
interchange.collections['vdW'].potentials=water_interchange.collections['vdW'].potentials
interchange.collections['Electrostatics'].potentials=water_interchange.collections['Electrostatics'].potentials

interchange.to_prmtop("water.prmtop")
interchange.to_inpcrd("water.crd")

Output

With time (BLUE) image

After deleting time (RED) image

GregorySchwing commented 10 months ago

VMD appears to read rst7 and crd files differently.

The time value messed up both, but if you remove the time column the crd files are still not read correctly. The rst7 file looks right.

Here is the exact same file with the extension changed from crd to rst7

image

Not sure this is an interchange problem, so you can close if you like

GregorySchwing commented 9 months ago

CRD files generated by Tleap and interchange tkdiff'ed. TIP3 TLEAP image

TIP3 Interchange image

mattwthompson commented 9 months ago

The time value isn't doing much useful in these files (in restart files they'd be important) but as far as I can tell they're part of the format. Here's what I'm going off of this documentation: https://ambermd.org/FileFormats.php#restart

I don't know what's supposed to be in an rst7 file, I think sometimes an .inpcrd is synonymous but it's opaque to me. If this is documented somewhere I'm open to adding a dedicated writer, but as I understand the ecosystem now I think just writing out .inpcrd files is all that Interchange needs to do