Moltemplate: Number of dihedrals in data file output

Sean · December 20, 2023, 2:33am

Hello,

I recently built a system with moltemplate, and I noticed a discrepancy between the number of dihedrals that the LAMMPS Data file says are present and the number of dihedrals that should be present based on my calculations, which was merely counting the number of defined dihedrals within each monomer (I have a two monomer copolymer), counting the number of defined dihedrals between each monomer (each polymer chain consists of 24 repeat units), and then summing them for every chain in the system.

Using this counting method, I’m getting that there are 300 more dihedrals defined in the system than there are according to the LAMMPS data file produced by moltemplate. I have not been able to figure out why there are more based on counting and, if some are missing in the as-built system, which ones.

Has anyone come across a problem like this? For the person so inclined, I have attached all of the .lt files for the moltemplate build as well as a summary of my counting. Any help is much appreciated.

Kind regards,
Sean

Moltemplate dihedrals.xlsx (10.4 KB)
AAPoly_175_10242023.zip (397.5 KB)

akohlmey · December 20, 2023, 11:26pm

When using the VMD TopoTools scripts that derive angles and dihedrals from the bond topology, I get the following numbers.

Original data file:

Info) Analyzing structure ...
Info)    Atoms: 9675
Info)    Bonds: 10550
Info)    Angles: 10125  Dihedrals: 21175  Impropers: 0  Cross-terms: 0
Info)    Bondtypes: 38  Angletypes: 7  Dihedraltypes: 12  Impropertypes: 0
Info)    Residues: 650
Info)    Waters: 0
Info)    Segments: 1
Info)    Fragments: 25   Protein: 0   Nucleic: 0

After regenerating angles and dihedrals:

Info) Analyzing structure ...
Info)    Atoms: 9675
Info)    Bonds: 10550
Info)    Angles: 16825  Dihedrals: 25200  Impropers: 0  Cross-terms: 0
Info)    Bondtypes: 38  Angletypes: 11  Dihedraltypes: 17  Impropertypes: 0
Info)    Residues: 650
Info)    Waters: 0
Info)    Segments: 1
Info)    Fragments: 25   Protein: 0   Nucleic: 0

I could only speculate on possible reasons but only @jewettaij as the author would know for certain.

Sean · December 21, 2023, 2:50am

Okay, thank you for checking, Axel.

hothello · December 21, 2023, 4:26pm

When a list of bonds, angles and dihedrals is explicitly set in the LT files, Moltemplate generates a molecular topology that is consistent with (and restricted to) the input terms. The discrepancy that you observe originates from missing terms in the monomers and, I suspect, in the polymer topology. Take the file HBA_RU_10132023.lt:

	Bonds	Angles	Dihedrals
LT	14	11	22
TopoTools	13	19	24

In the polymer this difference adds up. I have defined an oligomer to help debug your system:

test01.lt

# Original force filed.
import "OPLS_HBAHNA.lt"
import "HBA_RU_10132023.lt"
import "HNA_RU_10132023.lt"
import "EndgroupH.lt"
import "EndgroupOH.lt"

poly inherits OPLSAAmod {
 Endgrouper  = new EndgroupH.move(.5,0,-1)
 monomersHBA = new HBA_RU
 monomersHNA = new HNA_RU.move(0,0,5)
 Endgroupee  = new EndgroupOH.move(0,0,15.5)
 
 write('Data Bonds') {
   $bond:o3_EGh  @bond:o3EGh  $atom:Endgrouper/h1   $atom:monomersHBA/o3
   $bond:c7o8_1  @bond:c7o8   $atom:monomersHBA/c7  $atom:monomersHNA/o8
   $bond:c16_EGo @bond:c16EGo $atom:monomersHNA/c16 $atom:Endgroupee/o1
 }
 write("Data Angles"){
   $angle:c7_o8_c9_0   @angle:90_408_406 $atom:monomersHBA/c7   $atom:monomersHNA/o8  $atom:monomersHNA/c9
   $angle:c21_c7_o8_0  @angle:90_90_408  $atom:monomersHBA/c21  $atom:monomersHBA/c7  $atom:monomersHNA/o8
   $angle:c6_c7_o8_0   @angle:90_90_408  $atom:monomersHBA/c6   $atom:monomersHBA/c7  $atom:monomersHNA/o8
   $angle:c7_o8_c9_11  @angle:90_408_406 $atom:monomersHBA/c7   $atom:monomersHNA/o8  $atom:monomersHNA/c9
   $angle:c21_c7_o8_11 @angle:90_90_408  $atom:monomersHBA/c21  $atom:monomersHBA/c7  $atom:monomersHNA/o8
   $angle:c6_c7_o8_11  @angle:90_90_408  $atom:monomersHBA/c6   $atom:monomersHBA/c7  $atom:monomersHNA/o8
 }
 write("Data Dihedrals"){
    $dihedral:c7_o8_c9_o10_0  @dihedral:90_408_406_407  $atom:monomersHBA/c7  $atom:monomersHNA/o8  $atom:monomersHNA/c9  $atom:monomersHNA/o10
    $dihedral:c6_c7_o8_c9_0    @dihedral:90_90_408_406    $atom:monomersHBA/c6  $atom:monomersHBA/c7  $atom:monomersHNA/o8  $atom:monomersHNA/c9
    $dihedral:c5_c6_c7_o8_0    @dihedral:90_90_90_408    $atom:monomersHBA/c5  $atom:monomersHBA/c6  $atom:monomersHBA/c7  $atom:monomersHNA/o8
    $dihedral:c21_c7_o8_c9_0  @dihedral:90_90_408_406    $atom:monomersHBA/c21  $atom:monomersHBA/c7  $atom:monomersHNA/o8  $atom:monomersHNA/c9
    $dihedral:c22_c21_c7_o8_0  @dihedral:90_90_90_408    $atom:monomersHBA/c22  $atom:monomersHBA/c21  $atom:monomersHBA/c7  $atom:monomersHNA/o8
    $dihedral:c7_o8_c9_o10_11  @dihedral:90_408_406_407  $atom:monomersHBA/c7  $atom:monomersHNA/o8  $atom:monomersHNA/c9  $atom:monomersHNA/o10
    $dihedral:c6_c7_o8_c9_11  @dihedral:90_90_408_406    $atom:monomersHBA/c6  $atom:monomersHBA/c7  $atom:monomersHNA/o8  $atom:monomersHNA/c9
    $dihedral:c5_c6_c7_o8_11  @dihedral:90_90_90_408    $atom:monomersHBA/c5  $atom:monomersHBA/c6  $atom:monomersHBA/c7  $atom:monomersHNA/o8
    $dihedral:c21_c7_o8_c9_11  @dihedral:90_90_408_406    $atom:monomersHBA/c21  $atom:monomersHBA/c7  $atom:monomersHNA/o8  $atom:monomersHNA/c9
    $dihedral:c22_c21_c7_o8_11  @dihedral:90_90_90_408    $atom:monomersHBA/c22  $atom:monomersHBA/c21  $atom:monomersHBA/c7  $atom:monomersHNA/o8
 }
}

polymer = new poly [1]

write_once("Data Boundary") {
     0 20  xlo xhi
     0 20  ylo yhi
   -10 20 zlo zhi
}

With your original force field definition, this is what VMD TopoTools see:

Info) Analyzing structure ...
Info)    Atoms: 35
Info)    Bonds: 38
Info)    Angles: 31  Dihedrals: 66  Impropers: 0  Cross-terms: 0
Info)    Bondtypes: 36  Angletypes: 7  Dihedraltypes: 12  Impropertypes: 0

While after regenerating angles and dihedrals:

Info) Analyzing structure ...
Info)    Atoms: 35
Info)    Bonds: 37
Info)    Angles: 57  Dihedrals: 84  Impropers: 0  Cross-terms: 0
Info)    Bondtypes: 8  Angletypes: 11  Dihedraltypes: 17  Impropertypes: 0
Info)    Residues: 4

Note that in your forcefield 7 unique angle types are present, vs 11 found by TopoTools. Same story for the dihedrals, plus one extra bond defined in the monomer HBA. I have re-defined the topology of the monomers using the auto-generation of angle and dihedral terms. I have also cleaned up duplicate bond, angle, and dihedral coefficients, and defined 4 angle types that were not present in your force field. After the following, the angle count is correct:

Info) Analyzing structure ...
Info)    Atoms: 35
Info)    Bonds: 37
Info)    Angles: 57  Dihedrals: 68  Impropers: 0  Cross-terms: 0
Info)    Bondtypes: 7  Angletypes: 11  Dihedraltypes: 7  Impropertypes: 0
Info)    Residues: 4

I have no reason to believe the dihedral count is buggy in Moltemplate. I leave to you to work out the 10 missing dihedral types, that should be added to the file OPLS_HBAHNA_rules.lt. This is the oligomer with automatically generated angle and dihedral terms:

test03.lt

# Let Moltemplate guess the list of angles and dihedrals.

import "OPLS_HBAHNA_rules.lt"
import "HBA_RU_rules.lt"
import "HNA_RU_rules.lt"
import "EndgroupH.lt"
import "EndgroupOH.lt"

poly inherits OPLSAAmod {
	Endgrouper  = new EndgroupH.move(.5,0,-1)
	monomersHBA = new HBA_RU
	monomersHNA = new HNA_RU.move(0,0,5)
	Endgroupee  = new EndgroupOH.move(0,0,15.5)

write('Data Bond List') {
		$bond:o3_EGh  $atom:Endgrouper/h1   $atom:monomersHBA/o3
		$bond:c7o8_1  $atom:monomersHBA/c7  $atom:monomersHNA/o8
		$bond:c16_EGo $atom:monomersHNA/c16 $atom:Endgroupee/o1
	}
}

polymer = new poly [1]

write_once("Data Boundary") {
	   0 20  xlo xhi
	   0 20  ylo yhi
	 -10 20 zlo zhi
}

and the modified files:
lt_with_rules.tgz (4.9 KB)

@akohlmey I have used the following commands to recompute the angles and dihedrals:

topo readlammpsdata test01.data full
topo clearangles
topo cleardihedrals
topo guessangles
topo guessdihedrals
topo writelammpsdata test01.data2

Sean · December 22, 2023, 7:10pm

Thank you for your time on this - this is very helpful. I’m digesting this and debugging with the files you provided.

hothello · December 22, 2023, 11:06pm

Your’re welcome. You can already define the dihedrals rules from the existing coefficients, e.g.:

write_once("In Settings"){
# dihedral_coeff  @dihedral:90_90_406_408B  -0.123778   16.078802    0.251615    -1.616329     # Not used (same rule)
# dihedral_coeff  @dihedral:90_90_406_408N  -0.07425369 16.60391196  0.393870025 -1.695749715  # Not used (same rule)
  dihedral_coeff  @dihedral:90_90_406_408    0          2.1          0            0
  dihedral_coeff  @dihedral:407_406_90_90    0          2.1          0            0
  dihedral_coeff  @dihedral:406_90_90_90     0          7.25         0            0
  dihedral_coeff  @dihedral:90_90_90_90      0          7.25         0            0
  dihedral_coeff  @dihedral:90_90_90_91      0          7.25         0            0
  dihedral_coeff  @dihedral:406_90_90_91     0          7.25         0            0
  dihedral_coeff  @dihedral:408_90_90_91     0          7.25         0            0
  dihedral_coeff  @dihedral:90_90_90_408     0          7.25         0            0
  dihedral_coeff  @dihedral:90_90_408_406    0          2.5          0            0
# dihedral_coeff  @dihedral:90_90_406_407    0          2.1          0            0 # Duplicate of dihedral:407_406_90_90
# dihedral_coeff  @dihedral:90_90_90_406     0          7.25         0            0 # Duplicate of dihedral:406_90_90_90
  dihedral_coeff  @dihedral:90_408_406_407   0          5.0          0            0
  dihedral_coeff  @dihedral:91_90_90_91      0          7.25         0            0
}

From which the following 11 rules are derived:

  write_once("Data Dihedrals By Type") {
    @dihedral:406_90_90_90   @atom:406 @atom:90  @atom:90  @atom:90
    @dihedral:406_90_90_91   @atom:406 @atom:90  @atom:90  @atom:91
    @dihedral:407_406_90_90  @atom:407 @atom:406 @atom:90  @atom:90
    @dihedral:408_90_90_91   @atom:408 @atom:90  @atom:90  @atom:91
    @dihedral:90_408_406_407 @atom:90  @atom:408 @atom:406 @atom:407
    @dihedral:90_90_406_408  @atom:90  @atom:90  @atom:406 @atom:408
    @dihedral:90_90_408_406  @atom:90  @atom:90  @atom:408 @atom:406
    @dihedral:90_90_90_408   @atom:90  @atom:90  @atom:90  @atom:408
    @dihedral:90_90_90_90    @atom:90  @atom:90  @atom:90  @atom:90
    @dihedral:90_90_90_91    @atom:90  @atom:90  @atom:90  @atom:91
    @dihedral:91_90_90_91    @atom:91  @atom:90  @atom:90  @atom:91
  }

This leaves out 6 terms for which you need to define coefficients and rules. Hint: I have used the data file generated by TopoTools to see the list of dihedrals.

@atom:90  @atom:90  @atom:90  @atom:109
@atom:90  @atom:90  @atom:109 @atom:110
@atom:90  @atom:406 @atom:408 @atom:90
@atom:90  @atom:406 @atom:408 @atom:212
@atom:91  @atom:90  @atom:90  @atom:109
@atom:407 @atom:406 @atom:408 @atom:212

I hope this helps.

PS don’t forget to cite LAMMPS and Moltemplate in the resulting publication

Sean · December 23, 2023, 1:08am

I will certainly remember to cite both LAMMPS and Moltemplate! Couldn’t have done it without them

I wanted to address some of the items you pointed out:

The 4 unique missing angle types are end group angles or angles that terminate in aromatic hydrogens, which I left unspecified on purpose but, you are correct, are missing. I naively assumed they would treated with some default values, but I actually am not sure how lammps treats undefined angles, etc.
In the OPLS_HBAHNA_rules file, I went through just now and took out the ones that were indeed redundant and compared with my OPLS_HBAHNA file after doing the same. The only dihedral that I don’t have is the one terminating in an aromatic hydrogen on each end, which I think I knew at one time that I had left out on purpose–but again you are correct that it was missing. I also did not bother–perhaps wrongly–to define dihedrals involving the end groups. That accounts for 5 of the 6 dihedrals for which I do not have types defined. The last is what I’ve termed for myself a “spanning dihedral” since it spans from one monomer to another, and may (?) not be necessary if the other dihedrals controlling rotation about those bonds are present. But, again, it is missing as you said.
And yes, you are correct again, the way I defined bond types is redundant and unnecessary. I don’t recall the initial reason for doing it that way, but it is indeed a rookie mistake.
There is also a duplicate bond, line 47 of my HBA_RU_10132023.lt file, which accounts for the 38 bonds versus the 37 in your corrected files.
The separate 90_90_406_408B and 90_90_406_408N dihedrals I specify specifically rather than using a single 90_90_406_408 because I did quantum mechanical conformational calculations on those different dihedrals, and I am using fits to the potential energy curves from those calculations rather than the coefficients in the OPLS distribution.

So, I agree there were definitely errors in my files, and I am grateful for you pointing them out.

What I’m still trying to understand is, even with my duplicate dihedral types, if I asked moltemplate to create 19,250 actual dihedral instances and it only created 19,050–where are the 200 dihedrals I asked to be specified? And which dihedrals did it choose not to specify? I get the numbers 19,050 and 19,250 from the lammps data file generated by my original moltemplate input files and from my calculations based on the number of dihedrals within and between monomers, respectively. It could possibly be that my counting is wrong, but as I’ve thought about it I don’t see where I might have gotten off.

If I am missing something obvious, I apologize! Just trying to figure out where I’ve gone wrong. Your help is so much appreciated!

hothello · December 23, 2023, 9:11am

Regarding your points:

LAMMPS only does what you tell the program to do. If some terms are missing in the molecular topology (e.g. angle and dihedral terms in the data file), then there won’t be any potential restricting the motion of those atoms. What a potential does is to correlate the motion of groups of atoms, but LAMMPS doesn’t care whether those terms are specified or not.
You wonder why the dihedral counts was lower in your model, and I explained why:

Moltemplate generates a molecular topology that is consistent with (and restricted to) the input terms.

By using the write("Data Angles"){} and write("Data Dihedrals") {} blocks, you restrict Moltemplate to use only the specified bonds.
When you ask Moltemplate to use rules for the creation of angle and dihedrals, you won’t be able to distinguish between terms that involves the same atom types (e.g. 1-1-6-7) in different parts of a molecule. To be able to make this distinctions, you need to set all those terms explicitly and to assign different coefficient-type variables @- to each term.
Sean:

if I asked moltemplate to create 19,250 actual dihedral instances and it only created 19,050–where are the 200 dihedrals I asked to be specified?

This is not correct. You asked Moltemplate to create a polymer using:
- 12 HBA monomers, each with 22 dihedrals.
- 12 HNA monomers, each with 40 dihedrals
- 115 extra dihedral terms in the polymer.
  The math is: 12*22+12*40+115=859. Note that by default, Moltemplate will remove any duplicate terms unless you specify the options -overlay-bonds -overlay-angles -overlay-dihedrals, so the number you get is 847 because the residue HNA has one duplicate dihedral term between atoms 13 4 2 1. I have introduced the shortcut option -overlay-all to remove checks on all bonded terms, but made a syntax error in BASH. This is now corrected and I have opened a pull request.

The bottom line is:

moltemplate.sh -overlay-all sample.lt → will keep the list of bonded terms specified in the LT files. It allows to overlap multiple terms for the same bond/angle/dihedral/improper term.

moltemplate.sh sample.lt → will remove any duplicate bonded term, including angles, dihedrals, and impropers.

In both cases, the number bonded terms computed by Moltemplate won’t change if rules are used to create the molecular topology starting from the list of bonds.

I hope this is clear. It was useful to check the code again, cheers!

Sean · December 23, 2023, 1:18pm

Yes, this definitely helps! Thanks again!

Sean · January 2, 2024, 8:21pm

Still sorting through this - as a follow-up, I noticed in the output_ttree directory of a different build (that used the same monomer and forcefield input files), in the dihedral template file (attached), that half the number of a certain dihedral I defined is missing.

The dihedral is

$dihedral:c20c11c9o8	@dihedral:90_90_406_408N

It is defined in my monomer file HNA_RU_10132023.lt. It is right next to another 408N dihedral that does show up in the template file:

$dihedral:c12c11c9o8	@dihedral:90_90_406_408N

These missing dihedrals account for the difference in dihedrals between what’s given in the lammps data file and the number that should be there based on what was specified in the input files.

Is it normal for the ttree template files to not contain all of the dihedrals? Not trying to beat a dead horse, but I would like to understand the difference between what’s in my input files and what’s in the ttree dihedral template file. I’ve included all the files for re-creating my system (a random copolymer of 17 HBA, 7 HNA units per chain, 25 chains total).

Thanks!

Data Dihedrals.template (4.2 MB)
EndgroupH.lt (333 Bytes)
EndgroupOH.lt (512 Bytes)
HNA_RU_10132023.lt (8.1 KB)
HBA_RU_10132023.lt (4.9 KB)
OPLS_HBAHNA.lt (4.1 KB)
polybuild1_10132023_test.lt (30.9 KB)

Please note: These are the original files that contain redundant lines. Shouldn’t affect this specific missing dihedral question.

PS. Thought I should add the following: The 90_90_406_408N dihedrals only show up within HNA monomers, defined in the HNA_RU lt file. There are 2 of these 408N dihedrals defined per HNA monomer. There are 7 monomers per chain and 25 chains total = 2 x 7 x 25 = 350 total 408N dihedrals that should be present. In the dihedrals template file in the ttree directory, it only has 1 x 7 x 25 = 175 of these 408N dihedrals.

hothello · January 3, 2024, 10:17am

I am not sure I have the energy to debug your files further. As I have pointed out, Moltemplate will check and remove redundant terms unless the options -overlay-bonds -overlay-angles -overlay-dihedrals -overlay-impropers are used. With those options, you should get the exact number of terms specified in the monomer and polymer LT files. If you want all the terms specified, then you should use the rule-based generation of the molecular topology and specify the coefficients between atom types to ensure that every combination is covered. Downside: this method won’t allow to distinguish between dihedral terms involving the same atom types but located in different portions of a large molecule, which may benefit from a custom potential.

Sean · January 3, 2024, 7:47pm

Understood. Thank you for your time!