NAN when running GAMD with hydrogen mass repartitioned

I am getting a NaN error when trying to run gaMD on a PDB that has been prepped with heavy hydrogens.

Below are the XML scripts, the output, and the openMM prep script:
```xml
<?xml version="1.0" ?>
<gamd>
    <temperature>300</temperature> 

    <system>
        <nonbonded-method>PME</nonbonded-method>
        <nonbonded-cutoff>1.0</nonbonded-cutoff> 
        <constraints>AllBonds</constraints>
    </system>

    <barostat>
        <pressure>1.0</pressure> 
        <frequency>25</frequency>
    </barostat>
    
    <run-minimization>True</run-minimization>

    <integrator>
        <algorithm>langevin</algorithm>
        <boost-type>lower-dual</boost-type>
        <sigma0>
            <primary>6.0</primary> 
            <secondary>6.0</secondary> 
        </sigma0>
        <random-seed>0</random-seed>
        <dt>0.004</dt> 
        <friction-coefficient>1.0</friction-coefficient> 
        <number-of-steps> 
            <conventional-md-prep>2000</conventional-md-prep>
            <conventional-md>10000</conventional-md>
            <gamd-equilibration-prep>2000</gamd-equilibration-prep>
            <gamd-equilibration>20000</gamd-equilibration>
            <gamd-production>5000000</gamd-production> 
            <averaging-window-interval>50</averaging-window-interval>
        </number-of-steps>
    </integrator>

    <input-files>
      <forcefield>
        <coordinates>openmm_out.pdb</coordinates>
        <forcefields>
          <native>
            <file>amber99sb.xml</file>
            <file>tip3p.xml</file>
          </native>
        </forcefields>
      </forcefield>
      
    </input-files>

    <outputs>
        <directory>upper-dual/</directory>
        <overwrite-output>True</overwrite-output>

        <reporting>
            <energy>
                <interval>2000</interval>
            </energy>
            <coordinates>
                <file-type>DCD</file-type>
            </coordinates>
            <statistics>
                <interval>2000</interval>
            </statistics>
        </reporting>
    </outputs>
</gamd>
```

```
ERROR ~ Error executing process > 'run_gamd (1)'

Caused by:
  Process `run_gamd (1)` terminated with an error exit status (2)


Command executed:

  gamdRunner xml /home/delalamo/MD_nextflow/assets/gamd/gamd_test.xml
  python3 /home/delalamo/MD_nextflow/assets/gamd/process_traj.py -i openmm_out.pdb -t gamd_out/output.dcd -o out

Command exit status:
  2

Command output:
  windowCount:  0.0
  stage:  5.0
  stageOneIfValueIsZeroOrNegative:  31680864036.0
  stageTwoIfValueIsZeroOrNegative:  29910755042.0
  stageThreeIfValueIsZeroOrNegative:  28220817042.0
  stageFourIfValueIsZeroOrNegative:  24880955042.0
  stageFiveIfValueIsZeroOrNegative:  -275786936958.0
  Vmax_Dihedral:  5894.478369001709
  Vmax_Total:  -319653.50765317533
  Vmin_Dihedral:  5187.2246249725085
  Vmin_Total:  -330301.7906581727
  Vavg_Dihedral:  5766.738878097525
  Vavg_Total:  -320360.8485051052
  oldVavg_Dihedral:  0.0
  oldVavg_Total:  0.0
  sigmaV_Dihedral:  35.57305478270273
  sigmaV_Total:  176.82732694887434
  M2_Dihedral:  0.0
  M2_Total:  0.0
  wVavg_Dihedral:  0.0
  wVavg_Total:  0.0
  k_Dihedral:  0.001413919697763682
  k_Total:  9.39118541018011e-05
  k0prime_Dihedral:  3.907255752791146
  k0prime_Total:  2.1371961753692177
  k0doubleprime_Dihedral:  0.0
  k0doubleprime_Total:  0.0
  k0doubleprime_window_Dihedral:  0.0
  k0doubleprime_window_Total:  0.0
  boosted_energy_Dihedral:  nan
  boosted_energy_Total:  nan
  check_boost_Dihedral:  1.0
  check_boost_Total:  1.0
  threshold_energy_Dihedral:  5894.478369001709
  threshold_energy_Total:  -319653.50765317533
  thermal_energy:  2.494338785445972
  collision_rate:  1.0
  vscale:  0.9960079893439915
  fscale:  0.003992010656008516
  noisescale:  0.14097909017884955
  StartingPotentialEnergy_Dihedral:  nan
  StartingPotentialEnergy_Total:  nan
  ForceScalingFactor_Dihedral:  nan
  ForceScalingFactor_Total:  nan
  BoostPotential_Dihedral:  nan
  BoostPotential_Total:  nan
  k0_Dihedral:  1.0
  k0_Total:  1.0
  sigma0_Total:  25.104
  sigma0_Dihedral:  25.104

Command error:
  windowCount:  0.0
  stage:  5.0
  stageOneIfValueIsZeroOrNegative:  31680864036.0
  stageTwoIfValueIsZeroOrNegative:  29910755042.0
  stageThreeIfValueIsZeroOrNegative:  28220817042.0
  stageFourIfValueIsZeroOrNegative:  24880955042.0
  stageFiveIfValueIsZeroOrNegative:  -275786936958.0
  Vmax_Dihedral:  5894.478369001709
  Vmax_Total:  -319653.50765317533
  Vmin_Dihedral:  5187.2246249725085
  Vmin_Total:  -330301.7906581727
  Vavg_Dihedral:  5766.738878097525
  Vavg_Total:  -320360.8485051052
  oldVavg_Dihedral:  0.0
  oldVavg_Total:  0.0
  sigmaV_Dihedral:  35.57305478270273
  sigmaV_Total:  176.82732694887434
  M2_Dihedral:  0.0
  M2_Total:  0.0
  wVavg_Dihedral:  0.0
  wVavg_Total:  0.0
  k_Dihedral:  0.001413919697763682
  k_Total:  9.39118541018011e-05
  k0prime_Dihedral:  3.907255752791146
  k0prime_Total:  2.1371961753692177
  k0doubleprime_Dihedral:  0.0
  k0doubleprime_Total:  0.0
  k0doubleprime_window_Dihedral:  0.0
  k0doubleprime_window_Total:  0.0
  boosted_energy_Dihedral:  nan
  boosted_energy_Total:  nan
  check_boost_Dihedral:  1.0
  check_boost_Total:  1.0
  threshold_energy_Dihedral:  5894.478369001709
  threshold_energy_Total:  -319653.50765317533
  thermal_energy:  2.494338785445972
  collision_rate:  1.0
  vscale:  0.9960079893439915
  fscale:  0.003992010656008516
  noisescale:  0.14097909017884955
  StartingPotentialEnergy_Dihedral:  nan
  StartingPotentialEnergy_Total:  nan
  ForceScalingFactor_Dihedral:  nan
  ForceScalingFactor_Total:  nan
  BoostPotential_Dihedral:  nan
  BoostPotential_Total:  nan
  k0_Dihedral:  1.0
  k0_Total:  1.0
  sigma0_Total:  25.104
  sigma0_Dihedral:  25.104
```

```python3
#!/usr/bin/env python3

import argparse
import logging
import os
import sys
from sys import stdout

import openmm
import parmed
from openmm import app
from pdbfixer import PDBFixer

logging.basicConfig(level=logging.DEBUG)

def parse_args():
    """
    Parse command line arguments
    """
    argparser = argparse.ArgumentParser()
    argparser.add_argument("-i", "--input_pdb", required=True, type=str)
    argparser.add_argument("-o", "--output_pdb", default="openmm_out.pdb", type=str)
    argparser.add_argument("-v", "--verbose", action="store_true")
    return argparser.parse_args()

def fix_structure(pdbfile_in: str, pdbfile_out: str):
    """
    Fixed PDB by adding loops and replacing nonstandard AAs with standard 
    equivalents (for example, selenomethionine to methionine).
    Does not add hydrogens; this is done later.
    """
    fixer = PDBFixer(filename=pdbfile_in)
    fixer.findMissingResidues()
    fixer.findNonstandardResidues()
    fixer.replaceNonstandardResidues()
    fixer.removeHeterogens(False)
    fixer.findMissingAtoms()
    fixer.addMissingAtoms()
    app.PDBFile.writeFile(
        fixer.topology,
        fixer.positions,
        open(pdbfile_out, 'w'))

def solvate_protein(
        pdbfile_in: str,
        pdbfile_out: str,
        gmx_prefix: str="SYSTEM",
        forcefield: str="amber14-all.xml",
        water_model: str="amber14/tip3pfb.xml",
        ion_strength: float=0.15,
        positive_ion: str="K+",
        negative_ion: str="Cl-",
        padding: float=1.0,
        temperature: float=300.0,
        time_step_fs: float=4.0
        ):
    """
    Solvates model and adds ions.
    """
    pdb = app.PDBFile(pdbfile_in)
    forcefield = app.ForceField(forcefield, water_model)
    modeller = app.Modeller(pdb.topology, pdb.positions)
    modeller.addHydrogens(forcefield, pH=7.0)

    modeller.addSolvent(
        forcefield,
        padding=padding*openmm.unit.nanometers,
        ionicStrength=ion_strength*openmm.unit.molar,
        positiveIon=positive_ion,
        negativeIon=negative_ion)
    
    system = forcefield.createSystem(
        modeller.topology, 
        nonbondedMethod=app.PME,
        nonbondedCutoff=1*openmm.unit.nanometer,
        hydrogenMass=4.0 * openmm.unit.amu,
        constraints=app.AllBonds,
        rigidWater=False) # required to convert to gromacs

    integrator = openmm.LangevinIntegrator(
        temperature*openmm.unit.kelvin,
        1/openmm.unit.picosecond,
        time_step_fs*openmm.unit.femtoseconds
    )

    simulation = app.Simulation(
        modeller.topology,
        system,
        integrator
    )
    simulation.context.setPositions(modeller.positions)
    simulation.minimizeEnergy()
    positions = simulation.context.getState(getPositions=True).getPositions()

    app.PDBFile.writeFile(
        simulation.topology,
        positions,
        open(os.path.basename(pdbfile_out), 'w'))

    pmd_structure = parmed.openmm.load_topology(
        simulation.topology,
        system=system,
        xyz=positions)

    pmd_structure.save(f"{gmx_prefix}.top", overwrite=True)
    pmd_structure.save(f"{gmx_prefix}.gro", overwrite=True)

def main():
    """
    Main FXN
    """
    args = parse_args()
    pdbfile, out_filename = args.input_pdb, args.output_pdb
    if args.verbose: logging.basicConfig(level=logging.DEBUG)
        
    temp_pdb_filename = f"fixed_{os.path.basename(pdbfile)}"
    fix_structure(pdbfile, temp_pdb_filename)
    simulation = solvate_protein(temp_pdb_filename, out_filename)

if __name__ == "__main__":
    main()
```

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

NAN when running GAMD with hydrogen mass repartitioned #32

Metadata

Assignees

Labels

Projects

Milestone

Relationships

Development

NAN when running GAMD with hydrogen mass repartitioned #32

Description

Metadata

Metadata

Assignees

Labels

Projects

Milestone

Relationships

Development

Issue actions