[AMBER] Problem with Parmed-generated chamber prmtop

[hainm]

From amber mailing list, by Marc van der Kamp.

Just a placeholder, will update archive link with attached file later.

Original post

Dear all,

I have run into a curious problem when performing minimization with a FAD
co-factor, based on a Parmed-generated chamber prmtop.

Initially, a .psf file was generated with VMD's psfgen using CHARMM-style
parameters for FAD.
I confirmed that minimization with NAMD of FAD in vacuum, using the .psf
and .pdb file is ok.
(see min_namd.pdb attached).

Then I convert this psf and pdb with Parmed's chamber action (parmed
version 2.4.0 from AmberTools16):

chamber -top top_all36_prot.rtf -top top_all36_na.rtf -top
top_all36_cgenff.rtf -top jp4061522_si_002_mod.txt -param
par_all36_prot.prm -param par_all36_cgenff.prm -param par_all36_lipid.prm
-param jp4061522_si_003_mod.txt -str toppar_water_ions.str -psf
test_fad.psf -crd test_fad.pdb -box bounding

parmout test_fad.prmtop test_fad.rst

However, when I perform minimization with sander using the resulting
test_fad.prmtop and test_fad.rst files, something curious happens:
The central structure in the flavine ring distorts in a strange way - In
particular, C6A and C4A appear to exchange position, and both the
configuration around N5 and N10 is distorted.
(see min_test.pdb)

I've checked what parmed reports for the bonds of these atoms, and this is
all as expected.

Is this some strange bug in the writing of the prmtop?
I'm attaching an archive that contains the files required to run the above
Parmed-command.

Thank you in advance for looking into this!
Marc

[swails]

There's something wrong with improper torsions. This is what I saw when doing basic energy comparisons of the ParmEd-generated prmtop and the CHARMM psf with parameter files using OpenMM and sander API:

# OpenMM
In [28]: pmd.openmm.energy_decomposition_system(psf, psf_system)
Out[28]: 
[('HarmonicBondForce', 214.09383782289848),
 ('HarmonicAngleForce', 29.220837540853598),
 ('PeriodicTorsionForce', 53.14508795674756),
 ('HarmonicBondForce', 16.169664786221116),
 ('CustomTorsionForce', 0.10394766426839853),
 ('NonbondedForce', -196.16814733021687),
 ('CMMotionRemover', 0.0)]

In [29]: pmd.openmm.energy_decomposition_system(parm, parm_system)
Out[29]: 
[('HarmonicBondForce', 214.09383782289848),
 ('HarmonicAngleForce', 29.220837540853598),
 ('PeriodicTorsionForce', 53.14508795674756),
 ('HarmonicBondForce', 16.169664786221116),
 ('CustomTorsionForce', 0.10394766426839853),
 ('NonbondedForce', -196.16814733021687),
 ('CMMotionRemover', 0.0)]

# sander
In [37]: pmd.tools.energy(parm, igb=6).execute()
Bond          =          214.0940710     Angle         =           29.2204527
Dihedral      =           53.1450781     Urey-Bradley  =           16.1697296
Improper      =       274340.8989297     1-4 vdW       =           39.2153189     1-4 Elec.     =         -408.6796440
Lennard-Jones =           -4.7430512     Electrostatic =          178.0339330
TOTAL         =       274457.3548177

What you'll see here is that the OpenMM energy for both the Amber and CHARMM systems are identical. However, the sander energy has an absolutely absurd improper contribution. This smells strongly like a missing degree<->radian conversion.

So I went back and checked out the original chamber source code. On the one hand, I see this comment:

   write(unit,flag_format)"CHARMM_IMPROPER_PHASE"
   write(unit,'(a)')"%COMMENT  psi: degrees"
   write(unit,form_format)fmt
   write(unit,fmt)(phase_impr(i),i=1,nimprtypes)                                                                                                                                         

I read that as the value should be in degrees. But if you look at psfprm.F90, you see this:

     ! Rebuild improper list
     call allreal(ifound,pk_impr,phase_impr)                                                                                                                                             
     do n=1,nimprtypes
       pk_impr(n)=pk0(n,0)
       phase_impr(n)=pi*phase0(n,0)/180.d0
     enddo

I don't know... that looks to me an awful lot like the phase is being stored in radians, not degrees. This is a rookie mistake, and one that should never have made it this far. So why did it? Pretty simple -- 99.9% of all improper torsions have phase offsets of 0 -- the only angle with the same value measured in both degrees and radians. Marc happened to find the one improper torsion not like that.

This will be an easy fix, although quite annoying! I will have to fix the comment.