[AMBER] Fwd: problem with chamber & vmd generated psf (line 2349 of file psfprm.F90)

From: Sarah Witzke <witzke.sdu.dk>
Date: Thu, 7 Aug 2014 16:11:48 +0000

I just wanted to send this mail to the emailing list for the archives in case anyone should have trouble with missing bonds in TIP3 water from VMD.

Start på videresendt besked:

Fra: Michael Petersen <mip.sdu.dk<mailto:mip.sdu.dk>>
Emne: SV: [AMBER] problem with chamber & vmd generated psf (line 2349 of file psfprm.F90)
Dato: 7. aug. 2014 15.07.18 CEST
Til: Sarah Witzke <witzke.sdu.dk<mailto:witzke.sdu.dk>>, Marc van der Kamp <marcvanderkamp.gmail.com<mailto:marcvanderkamp.gmail.com>>

Hi Marc,

I think the problem with "fast 3-point water" arises because some of your water molecules lacks the HH bond usually employed to make it a rigid water molecule. My guess would be that you have added some water molecules using vmd (from e membrane patch possibly?).

If you find the solvate directory within the vmd tree (mine is /home/mip/vmd-1.9.2a39/plugins/noarch/tcl/solvate1.6) and look in wat.top you'll see that the line defining the H1 H2 bond is commented out. You need this bond to be defined if you want to run with pmemd/sander. namd does things slightly differently and implicitly adds the HH bond if it isn't defined so for most people using charmm it doesn't matter. But it does if you're using the amber MD engines.

You should also be able to see the HH bonds or lack of them in vmd if you load the psf file.

Hope this helps,
Michael

________________________________________
Fra: Sarah Witzke
Sendt: 7. august 2014 14:41
Til: Marc van der Kamp
Cc: Michael Petersen
Emne: Re: [AMBER] problem with chamber & vmd generated psf (line 2349 of file psfprm.F90)

Dear Mark,

One problem down - and then a new one arises - don't you just love research ;-)

Yes, I would prefer to change to ParmEd as well - also because Jason always keeps his code up to date and very user-friendly.

When the NBFIX terms came to our attention our solution was to edit the prmtop file (the one out of chamber) using (again) ParmEd and write a new .prmtop file. This is a bit of manual leg work that I wouldn't bother with when ParmEd's chamber is fully happy.
But just in case you want to know, I used the commands "changeLJPair .SOD @O13,014 3.16 -0.07502", "changeLJPair @SOD @O22,032 3.13 -0.07502" and "changeLJPair @SOD @CLA 3.731 -0.08388" where the atom names O13 and O14 correspond to atom types O2L and atom names O22 and O32 correspond to the type OBL.

It seems odd with the TIP3 water - especially since the psf file does indeed include the 3 bonds and also since sander does not complain until water number 10. It could not be related to a numbering problem in your input coordinates if your pdb file reaches residue numbers above 9999? Just a thought.
I have cc'ed Michael, he might have an idea about this.

Best regards,
Sarah

Den 07/08/2014 kl. 13.47 skrev Marc van der Kamp <marcvanderkamp.gmail.com<mailto:marcvanderkamp.gmail.com><mailto:marcvanderkamp.gmail.com>>:

Dear Sarah,

Thanks for these very useful instructions. I was able to generate a prmtop+inpcrd for my system with CGenFF atom types in this way - that is the benefit of having atom types as numbers instead of strings!
So this could be a general work-around (at least until chamber becomes more tolerant with long atom type names in XPLOR-style psf's).

I would still prefer to use parmed.py, as this will give the correct NBFIX's (as Jason indicated), but at least this is a start.

Now I 'just' have to solve the issue that sander reports:
Error: A residue defined as a "fast 3-point water"
       is not defined by a triangle of three bonds.
       Residue 527 contains 2 bonds.

Residue 527 would be 6th TIP3P water molecule in my system. I have confirmed that in the psf that I use as input for chamber, the 1st to 10th TIP3P molecule have the expected 3 bonds.

I'll probably need to work out how to check if this is also the case in the prmtop.
I may start a new thread on this, but I thought that Sarah or others may have come across similar errors?

Thanks!
Marc




On 6 August 2014 16:27, Sarah Witzke <witzke.sdu.dk<mailto:witzke.sdu.dk><mailto:witzke.sdu.dk>> wrote:
Dear Mark,

Oh yeah - I should be more careful when reading the mails ;-) Sorry, I get it now.
So, what I do for my systems to get the inpcrd and prmtop files:
1) Do whatever the modifications in VMD
2) Change the xplor-type psf file to a charmm type psf file with a tcl script like this:

xplor2charmm.psf:
****
package require psfgen
topology charmm_all36_280813.top
resetpsf
readpsf xplor.psf
writepsf charmm charmm.psf
exit
****

I run the script with "vmd -dispdev text -e xplor2charmm.tcl > xplor2charmm.out".

3) Insert a line of:
      0 0 !NUMLP NUMLPH

If you have no cmap term, then this line should be inserted at the end of the .psf file.
If you do have cmap terms, then this section should only be the second last section.
For example here is the last few lines of a psf file with cmap terms (!NCRTERM):

***
      1 0 !NGRP
      0 0 0

      0 0 !NUMLP NUMLPH


      3 !NCRTERM: cross-terms
     56 216 218 224 216 218 224 229
    239 241 243 261 241 243 261 263
    261 263 265 271 263 265 271 276
****

4) Run chamber (maybe ParmEd now)

I always have to manually insert the section with !NUMLP when using psf files from vmd in chamber.
Whether there is an easier way to do it, I don't know, but maybe this can help you too?

Best regards,
Sarah


Den 06/08/2014 kl. 14.01 skrev Marc van der Kamp <marcvanderkamp.gmail.com<mailto:marcvanderkamp.gmail.com><mailto:marcvanderkamp.gmail.com>>:

Dear Sarah,

Thanks for your effort in checking this.
To clarify:
I meant that I have not been able to generate prmtop+inpcrd for my 'real'
system of interest, for which I generated a psf with VMD's psfgen.
I have used chamber several times before (successfully) using psf+crd
generated by CHARMM. In my original post, I indicated that I was keen to
use chamber for psfgen generated psf's, now AmberTools14 update.7 included
"Better support for VMD psfgen files" (see
http://ambermd.org/bugfixesat.html).

I thought that this would mean that the X-PLOR style psf (with atom type
names instead of numbers, as you pointed out) would be accepted, but
perhaps I am wrong.
As far as I know, psfgen can only write X-PLOR style psf's.

Thanks,
Marc


On 6 August 2014 12:50, Sarah Witzke <witzke.sdu.dk<mailto:witzke.sdu.dk><mailto:witzke.sdu.dk>> wrote:

Hi Marc,

I am a bit confused - you said that you did not get a .inpcrd and a
.prmtop file yet, but further down in the mail you said that chamber did
successfully generate those files?
I have generated those files with chamber without problems - let me know
if you would like them.
Chamber does not like xplor type pst files due to the difference in the
5th column of !NATOM section - notice the atom names vs. a number:

non-xplor type:
    48 !NATOM
     1 A 1 MET N 72 -0.300000 14.0070


xplor type:
    48 !NATOM
     1 A 1 MET N NH3 -0.300000 14.0070


As for ParmEd I have not tried it - but I look forward to it :-) Thanks
Jason!

Best regards,
Sarah


Den 06/08/2014 kl. 09.55 skrev Marc van der Kamp <marcvanderkamp.gmail.com<mailto:marcvanderkamp.gmail.com><mailto:marcvanderkamp.gmail.com>
:

I'm attaching the tri-peptide test psf's, crd and pdb as described in my
previous message here (in test.tgz), in case anyone wants to test this.
Thanks,
Marc


On 6 August 2014 08:39, Marc van der Kamp <marcvanderkamp.gmail.com<mailto:marcvanderkamp.gmail.com><mailto:marcvanderkamp.gmail.com>>
wrote:

Dear Jason (and others),

I wasn't aware of this new functionality of ParmEd! Very impressed (and
great that there is the NBFIX option as well!) - many thanks for your
work,
Jason!

Now, I did encounter problems, and haven't been able to generate a
prmtop+inpcrd of my system yet; for a tri-peptide test-system, I was
only
successful with chamber (not parmed.py) and a CHARMM-style psf (written
by
CHARMM). See details below.

First, I believe the -cmap flag should not be used in parmed.py?
When I try (after cloning and installing parmed.py as per your
instructions - did this on my MacBook as away from work):
chamber -cmap -top top_all36_prot_ksi.rtf -param
par_all36_prot_cgenff_ksi.prm -psf ionized.psf -crd ionized.pdb
Action chamber failed
UnhandledArgumentWarning: -cmap

Without the -cmap flag but with specifying a -box, I get:
chamber -top top_all36_prot_ksi.rtf -param
par_all36_prot_cgenff_ksi.prm
-psf ionized.psf -crd ionized.pdb -box 72.18 83.12 69.77
Action chamber failed
InputError: Box must be 3 lengths or 3 lengths and 3 angles!

(and adding 90.0 90.0 90.0 for angles doesn't help)

How should I specify the 'boundingbox' option?

Then, I proceeded to test parmed.py's chamber command and the
$AMBERTOOLS/bin/chamber with a simple tri-peptide system. I generated 3
different psf files:
- test.psf: CHARMM-style psf with charmm (v36)
- test_xplor: X-PLOR-style psf as written by charmm (v36)
- dry.test.psf: X-PLOR style psf as written by VMD's psfgen

For the CHARMM-generated psf's, I used the corresponding test.crd file
written by charmm, for the psfgen-generated psf, I used the
corresponding
dry.test.pdb file written by VMD/psfgen. For the protein parameters and
topology, I used the most recent files available on MacKerell's website
(
from toppar_c36_dec13.tgz
<
http://mackerell.umaryland.edu/download.php?filename=CHARMM_ff_params_files/toppar_c36_dec13.tgz

).
Below are first the initial ATOM lines of the psf, and then the
'results'
from chamber and parmed.py's chamber command.

CHARMM-style psf:

   48 !NATOM
    1 A 1 MET N 72 -0.300000 14.0070 0
0.00000 -0.301140E-02
    2 A 1 MET HT1 32 0.330000 1.00800 0
0.00000 -0.301140E-02
    3 A 1 MET HT2 32 0.330000 1.00800 0
0.00000 -0.301140E-02

chamber:
chamber -cmap -top top_all36_prot.rtf -param par_all36_prot.prm -psf
test.psf -crd test.crd

(prmtop + inpcrd succesfully generated)


parmed.py:

chamber -top top_all36_prot.rtf -param par_all36_prot.prm -psf test.psf
-crd test.crd
Creating chamber topology file from PSF test.psf, RTF files
[top_all36_prot.rtf], and PAR files [par_all36_prot.prm] Coords from
test.crd. Using CMAP. GB Radius set mbondi.
Traceback (most recent call last):
File "/Users/marcvanderkamp/bin/parmed.py", line 161, in <module>
 parmed_commands.cmdloop()
File

"/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/cmd.py",
line 142, in cmdloop
 stop = self.onecmd(line)
File

"/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/cmd.py",
line 221, in onecmd
 return func(arg)
File "<string>", line 1, in <lambda>
File

"/Users/marcvanderkamp/Library/Python/2.7/lib/python/site-packages/ParmedTools/parmed_cmd.py",
line 141, in _normaldo
 action.execute()
File

"/Users/marcvanderkamp/Library/Python/2.7/lib/python/site-packages/ParmedTools/ParmedActions.py",
line 3563, in execute
 parm.load_coordinates(coords)
File

"/Users/marcvanderkamp/Library/Python/2.7/lib/python/site-packages/chemistry/amber/_amberparm.py",
line 940, in load_coordinates
 for i, atom in enumerate(self.atom_list):
AttributeError: 'ChamberParm' object has no attribute 'atom_list'



X-PLOR style psf written by CHARMM:

   48 !NATOM
    1 A 1 MET N NH3 -0.300000 14.0070 0
0.00000 -0.301140E-02
    2 A 1 MET HT1 HC 0.330000 1.00800 0
0.00000 -0.301140E-02
    3 A 1 MET HT2 HC 0.330000 1.00800 0
0.00000 -0.301140E-02

chamber:
chamber -cmap -top top_all36_prot.rtf -param par_all36_prot.prm -psf
test.psf -crd test.crd

At line 2349 of file psfprm.F90
Fortran runtime error: Bad value during integer read

parmed.py:
chamber -top top_all36_prot.rtf -param par_all36_prot.prm -psf
test_xplor.psf -crd test.crd
Creating chamber topology file from PSF test_xplor.psf, RTF files
[top_all36_prot.rtf], and PAR files [par_all36_prot.prm] Coords from
test.crd. Using CMAP. GB Radius set mbondi.
Traceback (most recent call last):
File "/Users/marcvanderkamp/bin/parmed.py", line 161, in <module>
 parmed_commands.cmdloop()
File

"/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/cmd.py",
line 142, in cmdloop
 stop = self.onecmd(line)
File

"/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/cmd.py",
line 221, in onecmd
 return func(arg)
File "<string>", line 1, in <lambda>
File

"/Users/marcvanderkamp/Library/Python/2.7/lib/python/site-packages/ParmedTools/parmed_cmd.py",
line 141, in _normaldo
 action.execute()
File

"/Users/marcvanderkamp/Library/Python/2.7/lib/python/site-packages/ParmedTools/ParmedActions.py",
line 3563, in execute
 parm.load_coordinates(coords)
File

"/Users/marcvanderkamp/Library/Python/2.7/lib/python/site-packages/chemistry/amber/_amberparm.py",
line 940, in load_coordinates
 for i, atom in enumerate(self.atom_list):
AttributeError: 'ChamberParm' object has no attribute 'atom_list'


X-PLOR-style psf written by VMD's psfgen:

   48 !NATOM
    1 PROA 1 MET N NH3 -0.300000 14.0070 0
    2 PROA 1 MET HT1 HC 0.330000 1.0080 0
    3 PROA 1 MET HT2 HC 0.330000 1.0080 0

chamber:
chamber -cmap -top top_all36_prot.rtf -param par_all36_prot.prm -psf
dry.test.psf -crd dry.test.pdb

At line 2349 of file psfprm.F90
Fortran runtime error: Bad value during integer read


parmed.py:
chamber -top top_all36_prot.rtf -param par_all36_prot.prm -psf
dry.test.psf -crd dry.test.pdb
Creating chamber topology file from PSF dry.test.psf, RTF files
[top_all36_prot.rtf], and PAR files [par_all36_prot.prm] Coords from
dry.test.pdb. Using CMAP. GB Radius set mbondi.
Action chamber failed
ChamberError: Problem reading CHARMM PSF


The parmed.py errors for the psf's written by CHARMM may be due to my
installation?
The parmed.py error for the psf written by psfgen appear to be a
formatting issue?

Any hints/help greatly appreciated, as always.

--Marc


On 5 August 2014 20:30, Jason Swails <jason.swails.gmail.com<mailto:jason.swails.gmail.com>> wrote:


On Aug 5, 2014, at 12:27 PM, Marc van der Kamp <
marcvanderkamp.gmail.com<mailto:marcvanderkamp.gmail.com>>
wrote:

Hello,

Encouraged by the recent update of chamber related to better support
for
VMD psfgen generated psf files, I decided to try it out (after
updating
AmberTools14 with updates 1-7 and recompiling).

Unfortunately, when I ran:
$AMBERHOME/bin/chamber -cmap -top top_all36_prot_lig.rtf -param
par_all36_prot_cgenff_lig.prm -psf ionized.psf -crd ionized.crd -p
1ohp_rs1_charmm.prmtop -inpcrd 1ohp_rs1_charmm.inpcrd -box 72.180
83.115
69.771

(or variations on this, e.g. with -str to add the additional
parameters/topogies)
I get the following error:

At line 2349 of file psfprm.F90
Fortran runtime error: Bad value during integer read

I've run out of time to do more testing today (and no time tomorrow),
so I
would like to check if anyone has seen this or has an idea what it may
be
caused by.

lines 2341-2349 in psfprm.F90:

    read(psf_unit,fmt01) ii, & ! 1
         lsegid, & ! AAL
         iresid(i), & ! 1
         lresat(i), & ! ASN
         atom_label(i), & ! N
         iac(i), & ! 54
         cg(i), & ! -0.47000
         amass(i), & ! 14.0070
         imove(i) ! 0


Any insight welcome; many thanks in advance!

I think we would actually need the PSF file in order to determine how
the
assumption of the line layout breaks the assumptions made in chamber.

There is an alternative to chamber that you can try, however. See
http://github.com/swails/ParmEd for an updated version of the ParmEd
program that is bundled with AmberTools. ParmEd has a new action
"chamber"
that will create chamber-style topology files from PSF and CHARMM
coordinate, restart, or PDB files with a syntax that is almost
identical to
the chamber program itself. I've found that it is quite a bit more
flexible in terms of what PSF files it will successfully parse, and it
has
a number of additional advantages as well. Particularly:

1. If you use VMD to create a PSF and PDB file, it will pull the box
information from the PDB file so you don't have to specify it in the
CHAMBER command.
2. It accepts a new keyword "boundingbox" that will assign an
orthorhombic box that encloses all atom centers
3. It supports and implements NBFIX modifications defined in CHARMM
parameter or stream files (ion parameters in CHARMM36 all have
NBFIXes, to
my knowledge)
4. The topology files generated by ParmEd are notably smaller because
it
compresses degenerate parameter types when possible.

Obviously #3 is the most important improvement, since it's the
difference
between "correct" and "incorrect" implementations of the CHARMM force
field
in the eyes of many reviewers, whereas #1, 2, and 4 are more
convenience
options than anything else (#4 can have performance implications,
particularly with pmemd.cuda when NBFIX mods are present. Fewer atom
types
means that the LJ coefficient matrices are more likely to fit inside
the
device cache and thereby avoid costly latencies associated with what is
effectively a cache miss -- I'm not sure when or if this becomes
significant or noticeable).

To use it, download and install ParmEd using a command like:

python setup.py build
python setup.py install [--user] [--install-scripts <directory>]

Where --install-scripts allows you to specify where parmed.py and
xparmed.py are installed, and the --user command will install the
Python
modules in your home directory (but still in a place where they will be
found by default by the Python interpreter). Then you can use
parmed.py
something like:

$ parmed.py
chamber -cmap -top top_all36_prot_lig.rtf \
       -param par_all36_prot_cgenff_lig.prm \
       -psf ionized.psf -crd ionized.crd \
       -box 72.180 83.115 69.771
outparm 1ohp_rs1_charmm.prmtop 1ohp_rs1_charmm.inpcrd

The only major differences in usage here is that you can supply
multiple
parameter or topology files, using a separate "-param
<parameter_file>",
"-top <rtf_file>", or "-str <stream_file>" argument. You can also add
the
keyword "usechamber" to have ParmEd invoke chamber to create the prmtop
file (although it won't work if chamber can't parse the PSF file).
This
was added to help people validate ParmEd's converter.


Another thing ParmEd can do is write PSF files. So if ParmEd can read
the PSF file, you can write a quick script that will read the PSF file
in
and then print a new one that chamber can recognize. A quick python
script
like:

from chemistry.charmm.psf import CharmmPsfFile

psf = CharmmPsfFile('ionized.psf')
psf.write_psf('fixed.psf')


will write a PSF file that chamber can read.

Hope this helps,
Jason

--
Jason M. Swails
BioMaPS,
Rutgers University
Postdoctoral Researcher
_______________________________________________
AMBER mailing list
AMBER.ambermd.org<mailto:AMBER.ambermd.org>
http://lists.ambermd.org/mailman/listinfo/amber
<test.tgz>_______________________________________________
AMBER mailing list
AMBER.ambermd.org<mailto:AMBER.ambermd.org>
http://lists.ambermd.org/mailman/listinfo/amber
_______________________________________________
AMBER mailing list
AMBER.ambermd.org<mailto:AMBER.ambermd.org>
http://lists.ambermd.org/mailman/listinfo/amber
_______________________________________________
AMBER mailing list
AMBER.ambermd.org<mailto:AMBER.ambermd.org><mailto:AMBER.ambermd.org>
http://lists.ambermd.org/mailman/listinfo/amber
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Thu Aug 07 2014 - 09:30:03 PDT
Custom Search