Preparing PDB file

Okay folks,
I use a combination of topotools and charmm2lammps now.

This seems to work,
if I'm missing a point,
I would be glad about advice.
Feel free to rap on my knuckles then.

Best,
S.

Okay folks,
I use a combination of topotools and charmm2lammps now.

This seems to work,
if I'm missing a point,
I would be glad about advice.
Feel free to rap on my knuckles then.

​using CHARMM with LAMMPS is a mixed blessing.
the CHARMM implementation in LAMMPS goes back a long time and seems to be
mostly based on CHARMM-19 and CHARMM-22.
also, LAMMPS was a much smaller package at the time a several "clean" and
practical features were invented after the basic charmm styles were
implemented and thus the code is a bit on the hackish side in how the
CHARMM specific force field features are implemented, e.g. the scaled 1-4
exclusions for dihedrals and dihedral multiplicity.

the syntax and structure of CHARMM force field and topology files has
evolved since, and it appears that CHARMM now also uses a slightly
different switching function for non-bonded interactions. if you want to do
proteins​, you should be aware, that the (important) CMAP correction is not
(yet) available in LAMMPS (people are working on it) and it is also
implemented in a somewhat unusual way (the facility to enable this has been
added to LAMMPS quite a while ago).

the code in charmm2lammps.pl has not seen much maintenance since a long
time ago and while it seems to be working well enough for many people, it
is unfortunately written in a style that makes maintenance of the code
extremely difficult.

assuming that you are ok with not reproducing the output of the current
CHARMm (code) with current CHARMM (parameters) *exactly*, the best
strategy is to use either the NAMD/psfgen toolchain or CHARMM-GUI or CHARMm
and then feed the resulting psf/pdb to charmm2lammps.pl and hope for the
best.

a more ambitious approach would be to implement the second step into
topotools. this should not be too difficult, as Josh Vermass from the
Schulten group at UIUC has contributed code for outputting gromacs topology
files for CHARMM force field simulations. programming a variant that would
output a LAMMPS data file with the parameters include or better yet a
corresponding include file for LAMMPS containing the required XXX_coeff
statements should be a doable project.

axel.

Dear Axel,

thanks for your detailed response.
I'm certainly for now okay with that strategy.

However, I might come back to you at some
point when I've time implementing the second
step into topotools directly.

Best,
Stephan

would the community benefit then more from a rewrite of charmm2lammps.pl
or adding forcefield output for LAMMPS from within the topotool package?

Best,
S.

> the code in charmm2lammps.pl <http://charmm2lammps.pl> has not seen much
> maintenance since a long time ago and while it seems to be working well
> enough for many people, it is unfortunately written in a style that
> makes maintenance of the code extremely difficult.

would the community benefit then more from a rewrite of charmm2lammps.pl
or adding forcefield output for LAMMPS from within the topotool package?

​if you would want to go as far as rewriting charmm2lammps.pl, i'd rather
suggest to implement this directly into LAMMPS as a "read_charmm" command
that would then take a psf format topology file, a charmm format parameter
stream and a coordinate/box dimension containing file as input.​ it should
just as easy to set whatever is needed directly into internal data
structures and you'd save having to output the data file, as LAMMPS already
has a write_data command.

but before that, i would rethink how to best realize the various functional
forms for bonded and non-bonded interactions.

...and at that point, one might consider to implement an even more radical
change by replacing the current numbered atom/bond/angle/dihedral/improper
types with symbolic names, which would make interfacing to force fields
like CHARMM, Amber, OPLS and others much more convenient.

axel.

Dear Axel,

thanks for all the suggestions.
I will need to think about this and then talk
with you about the new format to read into LAMMPS
probably at some point...

However I have one following question to that:

I write a LAMMPS data file with *topo*, however,
in my resulting data file, there is only 1 dihedral
and 1 improper type. This seems incorrect, right?

(I used *topo* from VMD 1.9.2 - so probably version 1.5
and I feed into VMD a PSF/PDB file pair)

Best,
S.

Additional information:
My PSF file I build with VMD too -
it uses the top_all27_prot_lipid_na.inp
topology file for that purpose.

Hope this is correct.
Then *topo* is applied to yield the LAMMPS data file.

Best,
S.

Dear LAMMPS users,

I tried another PDB file, however I receive
again only one angle and dihedral type in my output.

I'm kind of astonished.
I describe in detail again my procedure:

1. Download PDB from pdb.org as a .pdb file.
2. Load into VMD 1.9.2
3. Generate PSF file within VMD
4. Close VMD and open the PSF file
5. Load coordinates (PDB) into the PSF loaded already in VMD
6. Write LAMMPS data by topo writelammpsdata pdb.data full
(this is topotools 1.5)

Check the output - everything looks okay, except that only
one angle and dihedral type is used throughout the file.
After adding the appropriate Coeffs the simulations can
be started however the one angle and dihedral type is wrong.

Is there any obvious mistake I make?

best,
S.

Dear LAMMPS users,

I tried another PDB file, however I receive
again only one angle and dihedral type in my output.

I'm kind of astonished.
I describe in detail again my procedure:

1. Download PDB from pdb.org as a .pdb file.
2. Load into VMD 1.9.2
3. Generate PSF file within VMD
4. Close VMD and open the PSF file
5. Load coordinates (PDB) into the PSF loaded already in VMD
6. Write LAMMPS data by topo writelammpsdata pdb.data full
(this is topotools 1.5)

Check the output - everything looks okay, except that only
one angle and dihedral type is used throughout the file.
After adding the appropriate Coeffs the simulations can
be started however the one angle and dihedral type is wrong.

Is there any obvious mistake I make?

​yes. i think i have explained this a few times on this list and on vmd-l.

a psf file does not contain any type information for topology data and
topotools cannot write what it doesn't know.

to explain in more detail, in the charmm workflow this is not needed, as
the bond/angle/dihedral/improper types are inferred from the atom type
names and matched against the parameter file, which has the various
n-tuples of atom types​ (or wildcards).
a similar effect can be achieved in topotools using the
retypebond/retypeangles/retypedihedrals/... commands, however without any
wildcard matching.
this is where topotools and charmm2lammps.pl fundamentally differ.
topotools requires the user to provide all information that it should be
written to the data file, but is not restricted to a specific force field
family. please also note that topotools doesn't know anything about scaled
1-4 LJ parameters for dihedrals, which are defined in the parameter file as
well.

​axel.​

Thanks Axel!

I'm sorry that I've overseen your posts in the respective lists, I
apologize.

However your explanation helps me a lot,
so I know what to do now.

Best,
S.

Dear LAMMPS users,

I need to make sure that I've understood you correctly,
then I can probably write some code to automatize my task.

Topotools provides me with Bond, Angle, Dihedral and
Improper Coeffs hints in the written LAMMPS data file.
So I need to populate the Coeffs with values...

I would proceed as follows:
Given e.g. a 2-tuple (CC-CC) for one of the Bond Coeffs
I would look into the corresponding CHARMM parameter
file and try to find the line with the given 2-tuple,
this gives me the parameters I need write into my LAMMPS
data file for this 2-tuple. The same can be done for
Angle, Dihedrals and Improper Coeffs.

Is this correct?

Best,
Stephan

Dear LAMMPS users,

I need to make sure that I've understood you correctly,
then I can probably write some code to automatize my task.

Topotools provides me with Bond, Angle, Dihedral and
Improper Coeffs hints in the written LAMMPS data file.
So I need to populate the Coeffs with values...

I would proceed as follows:
Given e.g. a 2-tuple (CC-CC) for one of the Bond Coeffs
I would look into the corresponding CHARMM parameter
file and try to find the line with the given 2-tuple,
this gives me the parameters I need write into my LAMMPS
data file for this 2-tuple. The same can be done for
Angle, Dihedrals and Improper Coeffs.

Is this correct?

not entirely. in particular, not if you want to accurately represent
the CHARMM force field. CHARMM has some idiosyncrasies (like many
force fields of it kind):
- there are two sets of non-bonded parameters (PairCoeffs). normal and
for 1-4 interactions (i.e. atoms at the end of a dihedral angle). most
of the time they are the same, but some are scaled.
- there are overrides to the mixing rules for i,j nonbonded parameters
from the "NBFIX" section in the parameter stream files (e.g. for ions)
- dihedral parameter lists often include wildcards (atom type X), so
matching is a bit more complicated. the same dihedral might match
multiple times with different multiplicity in the dihedral interaction
for more complex dihedral potentials
- in CHARMM impropers are defined from the topology files. topotools
just guesses them for structures that look like they might be
requiring them to stay planar.
- CHARMM uses numerically fitted corrections for sequences of two
dihedral interactions, call cross-terms. people have been working on
implementing this in LAMMPS for years, but no code has yet been made
public.
- current CHARMM uses a slightly different switching function to
smoothly switch off LJ interactions than what LAMMPS implements. i
recall people mentioning on this list that they have an implementation
of them. they may even by code in the mailing list archives, but so
far none of that has migrated into the official LAMMPS code.
- please note that CHARMM requires using a modified TIP3P water
potential with LJ interactions on the hydrogen atoms.

to make certain, i strongly suggest to study the various publications
describing the CHARMM force field, its base parameter sets and how to
apply them, as well as making tests for individual components where
you compare against a code that fully implements the desired version
of the CHARMM force field to determine whether the differences are
significant for the simulations of interest or not. please note that
the CHARMM community is quite thorough in checking whether their force
field was applied correctly and thus any publication using CHARMM
might come under close scrutiny.

axel.

Dear Axel,
dear LAMMPS list,

this helps me a lot,
after some literature
research I could automate.

However, I encounter for instance
bond coefficient hints from Topotools,
e.g.

# C-CPH2
or # C-HC

which do *not* correspond to any pair
in the BONDS section of any CHARMM parameter
file. I'm not sure how to handle this.
Can somebody give me some suggestion?

Best regards,
Stephan

Dear Axel,
dear LAMMPS list,

this helps me a lot,
after some literature
research I could automate.

However, I encounter for instance
bond coefficient hints from Topotools,
e.g.

C-CPH2

or # C-HC

which do not correspond to any pair
in the BONDS section of any CHARMM parameter
file. I’m not sure how to handle this.
Can somebody give me some suggestion?

Are these atom names or atom types?
Are you able to run your system in charmm or namd?

Axel

Dear Axel,

these names refer to *atom types*,
I can run the simulation also in LAMMPS,
if I assign dummy coefficients for them.

Best,
Stephan

Dear Axel,

these names refer to *atom types*,
I can run the simulation also in LAMMPS,
if I assign dummy coefficients for them.

​it would be more important to know, if you can run a simulation with that
setup using CHARMM or NAMD without errors.​ if you can, then there is
something that you are overlooking (not sure what it is). if not, then you
have unusual or unexpected bonds. it is difficult to say so without knowing
your entire protocol.

but my gut feeling feeling is that you have unusual (=euphemism for wrong)
bonds. if i remember correctly, atom type "C" in CHARMM is the carbonyl
carbon and HC an N-terminal hydrogen. so that looks like your structure
building from the pdb wasn't quite right. the same goes for the other bond.

axel.

Dear Axel,
dear LAMMPS users,

I've one further question concerning
the PairCoeffs in LAMMPS.
You said the following:

- there are two sets of non-bonded parameters (PairCoeffs). normal and
for 1-4 interactions (i.e. atoms at the end of a dihedral angle). most
of the time they are the same, but some are scaled.

I wrote out easily the PairCoeffs,
however the dihedral_coeff required a value between (1.0 and 0.0)
in the very last column for each Dihedral Coefficient (weighting).

Here was told, that it would be reasonable to
set it to 1.0 for each Dihedral Coefficient:
http://lammps.sandia.gov/threads/msg24567.html

However, I'm not sure if this is still the suggestion.
If not I will accomodate my code.

Best regards,
Stephan

Dear Axel,
dear LAMMPS users,

I've one further question concerning
the PairCoeffs in LAMMPS.
You said the following:

> - there are two sets of non-bonded parameters (PairCoeffs). normal and
> for 1-4 interactions (i.e. atoms at the end of a dihedral angle). most
> of the time they are the same, but some are scaled.

I wrote out easily the PairCoeffs,
however the dihedral_coeff required a value between (1.0 and 0.0)
in the very last column for each Dihedral Coefficient (weighting).

Here was told, that it would be reasonable to
set it to 1.0 for each Dihedral Coefficient:
http://lammps.sandia.gov/threads/msg24567.html

However, I'm not sure if this is still the suggestion.

​what this has to be, is determined by the CHARMM parameter sets that you
are using.

in most cases, it will be 1.0, for some it will be 0.5 and some it wi​ll be
0.0.
i vaguely remember from discussions with josh vermaas, who added a CHARMM
export to gromacs topology files to topotools, that this is related to
cyclic systems, with 0.5 for 6-membered rings (as there 1-4 would be
counted twice with the native charmm way of determining 1-4 pairs) and 0.0
for 5 and smaller.

​axel.​