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.