How to create a periodic crystal structure with lots of angle, dihedral and improper interactions efficiently

LAMMPS LAMMPS Beginners
1

Dear lammps users:

I’m a theoretical physics graduate student who’s new to the Lammps software.
I want to reproduce the results of a paper(DOI:10.1103/PhysRevB.109.054103) which utilized lammps for hysteresis simulation of a spin crossover material.

The method used in the article is to perform all atom molecular dynamics simulation on a 101010 [Fe(pyrazine)] [Ni(CN)4] supercell(the lattice of [Fe(pyrazine)] [Ni(CN)4] looks like this[https://imgur.com/a/molecule-simulation-4fBYjoy], and the supercell looks like this[https://imgur.com/a/vn2ZKic]), which has 16000 atoms in total. This method requires assigning force field constants to the two-body bond/pair potential, three-body angle potential, four-body dihedral angle potential etc, and the force field structure looks like this[https://imgur.com/a/WhxyINV]. I want to produce this 101010 supercell model with the complete bond, angle, dihedral and improper etc information attached for lammps to use. In the process of trying to accomplish this task, I found the most difficult part is to assign the exact angle, dihedral and improper potentials to the right site of the structure, because in lammps one can only create one angle/dihedral/improper bond a time by using create_bonds command.

From what I have learned so far, I can think of 3 strategies to accomplish this:

a) <1> use “lattice custom” command to create a unit cell[https://imgur.com/a/molecule-simulation-4fBYjoy] , <2> “region xxx block 0 10 0 10 0 10” to create the 101010 supercell, <3>“create_atoms 1 region xxx basis 1 1…” to fill the lattice with atoms, <4> create two body bonds , <5> “write_data” to produce the data file, <6> assign angle, dihedral and improper potentials to the model one by one by editing the data file. # this is the most time consuming method which I definitely want to avoid.

b)<1> use “lattice custom” command to create a 221 lattice[https://imgur.com/a/6aHKnlB], <2>“region xxx block 0 1 0 1 0 1” to create a region, <3> “create_atoms 1 region xxx basis 1 1…” to fill the lattice with atoms, <4>create two body bonds, <5> “write_data” to produce the data file, <6> assign angle, dihedral and improper potentials to the model one by one by editing the data file, <7> “read_data” to read data file in a new in.x.lmp file, then use “replicate 5 5 10” command to create a 101010 supercell, <8> “delete_atoms overlap” and some other commands to delete the overlap atoms, add the missing bonds etc to produce the desired final model. # this method can save a lot of work. But the problem is that I don’t know how well and complete the “replicate” command could build the necessary bond, angle, dihedral and improper interactions between atoms. Besides, I heard that software like moltemplate etc can also replicate the structures. Does anyone have any idea which “replicate” could meet my needs better? Another problem is that it seems to me the “molecule” command can also produce a structure to be replicated, then does the “molecule” command suit me better than “lattice—region—…” method in this case?

c) <1> use “lattice custom” command to create a unit cell[https://imgur.com/a/molecule-simulation-4fBYjoy] , <2> “region xxx block 0 10 0 10 0 10” to create the 101010 supercell, <3>“create_atoms 1 region xxx basis 1 1…” to fill the lattice with atoms, <4> “fix bond/create” command to create bonds and “atype”, “dtype” and “itype” keywords to create the corresponding angle, dihedral and improper interactions. # the lammps documentation hasn’t talked about from what kind of rules the a,d and itype keywords create the angle, dihedral and improper interaction based on a pair of atoms, e.g. let’s say a bond is created between atom no. 1(belong to atom type a) and 2(type b) by command “fix bond/create … dtype 2 …”, then how can I be certain which other two atoms atom 1 and 2 create a dihedral bond with and at what kind of order? These rules are important in this case because the angle, dihedral and improper interactions have certain kind of structure[https://imgur.com/a/WhxyINV].

From your experiences, which kind of method is the better one to create the needed model in my case?

Comments and advice are very much appreciated, thank you!

2

There is an alternative way to create a viable topology for molecular crystals with bonds crossing the periodic boundaries, discussed in this mailing list thread and based on the software Moltemplate. To summarise:

  1. Create a supercell by replicating the unit cell (CIF or PDB) with your favourite tool. Mine is GDIS, available in almost every Linux distribution.
  2. Open the supercell (exported as a PDB) in VMD and use the topotools to create the bonds based on the distance.
  3. Export the structure plus bond connectivity to a LAMMPS data file.
  4. Convert the data file to a LT file with ltemplify.py.
  5. Add the force field terms for angles, impropers, and dihedrals to the LT file. Please refer to the “Bonded interactions by type” section of the Moltemplate manual.
    The method suggested implies familiarity with Moltemplate. Start by running some of the examples to get started, especially the ones on polymers.
1 Like
3

Thank you for your advice! I’ll check it out.
May I ask one question related to 4 body improper/dihedral interaction? I want to be sure how many improper/dihedral interactions there are for a certain kind of atoms set-up, e.g. in this picture(https://imgur.com/a/JGqMk5l, which is a screenshot from the force field structure[https://imgur.com/a/WhxyINV] given by the article the result of which I want to reproduce) an improper interaction is defined for these 5 atoms(red is atom type a, and green is atom type b), then how many improper interactions are there in this picture? is it one(the explicit dark line 1234) or two(1234 and 1254) or 4(1234 and 1325 and 1254 and 1543)?

4

Thank you very much for your advice, I finally completed the model, and from the result(https://imgur.com/a/yuGMCFK) of a test equilibrium run “velocity all create 1.0 87287 \ fix 1 all nvt temp 1.0 1.0 1 \ run 20000”, so far so good.

Before I completed the model by moltemplate, I had tried to construct the model by the method 2 mentioned above. I added by hand all the angle, dihedral and improper interactions to a 222 supercell(https://imgur.com/W8xYBVv). Then I used the “replicate 5 5 5 bond/periodic” command to make the 101010 supercell(https://imgur.com/U1A16gZ). Although, “replicate” command can indeed construct all the necessary two body bonds in the 101010 supercell, it definitely cannot construct all the angle, dihedral and improper interactions across boundary of the different replicated parts, which can be seen from the result(https://imgur.com/yLewG3R) of a trial test equilibrium run “velocity all create 1.0 87287 \ fix 1 all nvt temp 1.0 1.0 1 \ run 20000”. This conclusion can also been drawn if one looks into the angle, dihedral and improper data from the data file of the supercell given by the replicate command.

Although I have constructed a working model using moltemplate, there is one thing bugging me, the number(8000) of improper interactions generated by the moltemplate are actually less than the number(10500) I generated using the 2nd method. I don’t know for now if I added more improper interactions than needed by hand to the 222 supercell. Maybe I could get the answer when I finish the production run.

5

Hi Daniel,

I am glad to see that you managed to generate a valid input structure using Moltemplate. However, there should not be any need to specify the angles manually. Moltemplate can already work out the molecular topology (angles, dihedrals, and impropers) from the connectivity map. Moltemplate also looks (and removes) duplicate interaction, which would explain the difference you observe. Use the -overlay-impropers command to override this check, if this is consistent with the force field you are using.

1 Like
7

Thank you very much for the reply. In the last several days, I have looked into the discrepancy of the number of improper interactions generated by two methods. As the reply I submitted 12 days ago suggested, I added 4 improper interactions to the structure shown in this picture(https://imgur.com/a/JGqMk5l , I was unable to get any confirmation about how many improper interactions I need to add to this structure from the author who writes the paper) in the “replicate” method. However due to my ignorance I originally hadn’t used the appropriate codes to generate the desired improper interactions in the moltemplate method. As this exchange(Improper Dihedrals - #2 by Andrew_Jewett) shows, the first atom type in the moltemplate command “@improper:type1 @atom:typeA @atom:typeB @atom:typeC @atom:typeD”, which is typeA, refers to the so called central atom in the improper interaction, I didn’t know about this rule at first. As a result, certain kinds of improper interactions were not generated at all in my first attempt, which explains the cause of the discrepancy.

After correcting my codes, the number of improper interactions generated by the moltemplate method exceeded the replicate method in my latest attempt. And indeed, it’s a tricky issue to generate the right number of improper interactions in this kind of structure(https://imgur.com/a/JGqMk5l), in this case I need to designate two types of green atoms, green1 and green2 and sort atom types properly in moltemplate command to solve the problem.

As a side note, I found out due to different rules to define the angle between two planes adopted by different improper styles in lammps, one may not produce the final data structure readily used by lammps in one go. For instance, the so called central atom type may appear in the 2nd place in the lammps data structure, so one needs to exchange the 1st and 2nd column of the data generated by moltemplate.

8

Hi Daniel,

sorry for the late reply. I ran some tests to determine the behaviour of Moltemplate w.r.t. the order of atoms in the Impropers section.

The LAMMPS manual reports the expected behaviour:

Note that the ordering of the 4 atoms in an improper quadruplet determines the definition of the improper angle used in the formula for each style.

This means that the improper_style command sets the convention for the atom order. I checked this assumption on the benzene molecule as defined in the GROMOS-ATB force field, as it uses improper_style harmonic to control the motion of hydrogen atoms. Here are the atom names:
image

For the harmonic style, LAMMPS expects the 4 atoms in an improper quadruplet to be ordered I,J,K,L with bond connectivity I-J, I-K, I-L. This is how the first improper angle is defined in Moltemplate:

write("Data Atoms"){
 $atom:H1 $mol:... @atom:HC    0.126000 -2.3698428182E+00  7.4863313010E-01  1.7047870000E-04 # 1
 $atom:C2 $mol:... @atom:CAro -0.126000 -1.3329822200E+00  4.2113874790E-01 -4.7658900000E-05 # 2
 $atom:C3 $mol:... @atom:CAro -0.126000 -1.0312462267E+00 -9.4427126360E-01  2.1249608000E-03 # 3
 #...
}
write("Data Impropers"){
 $improper:i1 @improper:g1 $atom:C2 $atom:C11 $atom:H1 $atom:C3
 #...
}

And this is the corresponding entry DATA file:

Impropers

1 1 2 11 1 3
2 1 3 2 4 5
#...

If I swap the order of atoms I,J in $improper:i1, the corresponding entry in the DATA file is also swapped.

In conclusion, you don’t need to swap the order of atoms in post-processing but do so in the LT file where the improper angles are defined.

9

Hi hothello,

Thank you very much for your reply.

From my experience, I have a different opinion regarding this question. Please allow me to detail what I mean.

For the harmonic style, LAMMPS expects the 4 atoms in an improper quadruplet to be ordered I,J,K,L with bond connectivity I-J, I-K, I-L.

In LAMMPS documentation section 8.8, it states that the rule to define a harmonic improper angle is “…the angle between the plane of I,J,K and the plane of J,K,L.(ordered I,J,K,L)…” and “…Note that defining 4 atoms to interact in this way, does not mean that bonds necessarily exist between I-J, J-K, or K-L, as they would in a linear dihedral. Normally, the bonds I-J, I-K, I-L would exist for an improper to be defined between the 4 atoms…”. So I think atom I isn’t necessarily the so-called central atom.

But the default setting of improper interaction for moltemplate when using a user-specified force field without defining explicitly the rule of designating the central atom in the improper interaction I-J-K-L(which is my case) is that the central atom must be atom I in the moltemplate command “@improper:type1 @atom:typeI @atom:typeJ @atom:typeK @atom:typeL”.

For example, one of the improper interactions in the paper(10.1103/PhysRevB.109.054103) is “Fe–N–N–N (180◦)” as shown in this screenshot(https://imgur.com/a/improper-JGqMk5l 1 is Fe, 2/3/4 are N), the angle of which is between plane 134 and plane 123. So the moltemplate command write_once(“Data Impropers By Type”) must read “@improper:typeFeNNN @atom:Fe @atom:N @atom:N @atom:N”, which would produce a corresponding preliminary data file “Impropers 1 1 1 2 3 4”. However, the No. sequence in the final data file for the LAMMPS production run should be “Impropers 1 1 2 1 3 4”, which means I need to exchange the first and second column of the improper data matrix of the preliminary data file.

10

For working with Moltemplate and LAMMPS the correct ordering of atoms is the ordering that the software is programmed to use.

In this case it is having the improper’s central atom as the second atom in the list of four. As the LAMMPS documentation states (and thank you for checking it!!), this is so the improper angle can be defined as the angle between planes IJK and JKL (for an atom list IJKL). In turn, this simplifies the programming considerably since this is the same mathematical definition for a “proper” dihedral angle.

Thus, you should also list the central atom as the second of four in Moltemplate.

This is completely separate to how people list the atoms in a scientific paper meant for human readers. Neither Moltemplate not LAMMPS are reading the paper’s text when they run. Just your inputs.

11

Hi, srtee

Thank you very much for your reply.

Thus, you should also list the central atom as the second of four in Moltemplate.

This is completely separate to how people list the atoms in a scientific paper meant for human readers. Neither Moltemplate not LAMMPS are reading the paper’s text when they run. Just your inputs.

Regarding the above two comments, I want to clarify more clearly what I said.

From my email exchange with the author of moltemplate-Dr. Andrew Jewett, he stated:

Authors of the first molecular dynamics software programs did not bother to write separate routines for calculating dihedral forces and improper forces. They just re-used the dihedral-force code to implement improper forces. (Incidentally, for these forces, it is more common for the 2nd or 3rd atom to be the central atom. But there is no general consensus. When adding new force fields to moltemplate, this causes significant headaches. For the current implementation of OPLS, the 1st atom is the central atom. But for AMBER GAFF2, the 3rd atom is the central atom. See details below…)

  • The OPLS force field uses these rules (as indicated here). It generates only one improper interaction for every 4 atoms.

  • AMBER GAFF2 currently uses these rules (as indicated here). The 3rd atom is the central atom. For this force-field, moltemplates generates 3 improper interactions for every 4 atoms.

You can specify exactly how you want the impropers to be generated and what is the central atom. You can add your own “FILENAME.py” file to the “nbody_alt_symmetry/” folder and refer to it in parenthesis in the "Data Impropers By Type (FILENAME)

And according to my experience, when I don’t define explicitly the rule of designating the central atom in the improper interaction I-J-K-L of the user-specified force field, moltemplate wouldn’t generate any improper interaction if I didn’t place the central atom type at the 1st place of the moltemplate command write_once(“Data Impropers By Type”) “@improper:type1 @atom:A @atom:B @atom:C @atom:D”.

For example, in my last reply, if the moltemplate command write_once(“Data Impropers By Type”) read “@improper:typeFeNNN @atom:N @atom:Fe @atom:N @atom:N”, then moltemplate wouldn’t generate any typeFeNNN improper interaction.

12

If that was your actual input command then the reason for failure is more basic: an improper interaction must be created between atoms $atom: ... rather than atomtypes @atom: ....

13

In Moltemplate, bonded terms can be defined By Type and thus require @atomtype variables.

The problem reported by @DannLee is that Moltemplate by default assumes the central atom to be listed first, when it searches for bond topologies matching the specified atom types. For your command to work, you need to use the following command:

write_once("Data Impropers By Type (cenJsortIKL)"){
 @improper:typeFeNNN @atom:N @atom:Fe @atom:N @atom:N”
}

As Andrew said, the file in parenthesis specifies an alternative rule for searching the improper angles. Let us know if this work!

1 Like
14

Hi hothello,

After some experiments, I found out I was wrong about the feasibility of the LAMMPS command “replicate x x x bond/periodic” to produce a model of this periodic structure with all the necessary bond, angle, dihedral and improper interactions attached.

In my Jan 8 reply, I said

Before I completed the model by moltemplate, I had tried to construct the model by the method 2 mentioned above. I added by hand all the angle, dihedral and improper interactions to a 222 supercell(https://imgur.com/W8xYBVv). Then I used the “replicate 5 5 5 bond/periodic” command to make the 101010 supercell(https://imgur.com/U1A16gZ). Although, “replicate” command can indeed construct all the necessary two body bonds in the 101010 supercell, it definitely cannot construct all the angle, dihedral and improper interactions across boundary of the different replicated parts, which can be seen from the result(https://imgur.com/yLewG3R) of a trial test equilibrium run “velocity all create 1.0 87287 \ fix 1 all nvt temp 1.0 1.0 1 \ run 20000”. This conclusion can also been drawn if one looks into the angle, dihedral and improper data from the data file of the supercell given by the replicate command.

The above quoted conclusion is wrong. The reason I failed the last time is that I hadn’t added the angle, dihedral and improper interactions across the periodic boundary. Actually, it’s not even necessary to build a 222 supercell as an intermediate step. A unit cell data file with all the appropriate angle, dihedral and improper interactions manually added would do the job. After using the “replicate 10 10 10 bond/periodic” command to replicate the unit cell and then “write_data” to produce the data file, a final working model shall be available for a production run.

All in all, thank you very much for the valuable advice, your help is much appreciated:)

15

That’s good, but it still leaves the problem of how you count the 3- and 4-body interactions. It would be useful for me to test Moltemplate on the unit cell by specifying the bond connectivity and the higher bonded interactions By Type, with the additional symmetry rule for the improper angles.