Discontinuous density based on read_restart command

LAMMPS LAMMPS Mailing List Mirror
1

Dear lammps users,

I run a npt simulation based on a restart file trying to fully relax epoxy resin system and to calculate a reliable density value. However, at the beginning of the simulation there is a big discontinuity of the calculated density, i.e., the density shows a sudden increase compared to the previous values and then go stable slowly. As I can see, the density should be around the same value because restart files stored basically all information and there should not be too much difference between the new calculation and the previous one. Some of my code is below:

The previous calculation:

fix 1 all npt temp 500 500 100 iso 1 1 1000
velocity all create 500 8531416
compute 1 all pressure thermo_temp
fix 11 all ave/time 1 20 100 c_1[1] c_1[2] c_1[3] c_1[4] c_1[5] c_1[6] file anneal.data

compute 10 all msd
fix 91 all ave/time 1 20 100 c_10[4] file msd_500k.data
compute 11 all chunk/atom bin/1d x lower 1 units reduced
fix 92 all ave/chunk 100 10 1000 11 density/mass ave running file densit_500k.data

run 50000
unfix 1
write_restart equi500k.restart
uncompute 10
uncompute 11
unfix 91
unfix 92

2

The obvious things to check is to whether all pair styles you use write their state to the restart file and/or whether your kspace settings are the same. If that all checks out, then something might be amiss. Please also update to the latest version of Lammps so that we know you haven’t missed any bug fixes.

3

The obvious things to check is to whether all pair styles you use write
their state to the restart file and/or whether your kspace settings are the
same. If that all checks out, then something might be amiss. Please also
update to the latest version of Lammps so that we know you haven't missed
any bug fixes.

​i would also check whether the position and size of the bins remains
unchanged from before to after the restart.

axel.​

4

Dear Axel and Stefan,

Many thanks for the reply. I believe there should be some changes from before and after the restart. Due to my shortage of experience, could you roughly suggest how to check these changes ? (For example, the pair style state, the bin size and postiion ) Is all the informaton in the log files ? I also tried to check the restart file but seems the format can not be checked directly.

Regards

Hao

5

The restart file is binary so that is difficult to parse yourself. You should check the log files instead. Be sure to copy the old one or specify the log file with the “-log” command line switch so that it is not overwritten by the newer run.

6

Dear all,

I have been recently simulating the tensile behavior of highly crosslinked epoxy resin systems, in which bond breaking is necessary to be considered. I used fix bond/break command since it is the simplest way. However, atom missing happened sometime and I believe it is because of large forces generated when bond/angle/dihedrals suddenly disappeared when bond broke. I am also considering using fix bond/react instead but there are some issues not understood. The questions are as follows:

  1. for fix bond/break, what is the most effective way to prevent bond missing ? I am quite sure the topology is correct, and all measures I can take now is to reduce timesteps, enlarge cutoff/skin distance, perform comm_modify,etc. However since the large force due to sudden bond disappearing is unavoidable, all the above measures seem not able to provide a 100% guarantee.

  2. for fix bond/react, can it be used during tensile? i.e., combined with fix deform.

  3. still for fix bond/react, I sort of doubt if it can be used in my simulation. My system includes 60000 atoms. There are three bond types selected as the candidates to break , and the three bond types are successive in a chain. For example, in -O-C-C-N-, the three candidates are O-C, C-C and C-N, respectively (which bond break depends on the actual bond stretch length during tensile). I assume there should be three reactions taken into account which invovle fix bond/react three times in an input file. So a pre-reacted molecule template,a reacted template and a map file should be prepared for each of the candidate bond. However, assume one C-C bond is broken due to tensile (so the angle/dihedrals associated with this C-C is also deleted), it will make the template files for the O-C and C-N different. So how I can correctly prepare the template files in the first place?

Thanks for the patience. Just want to know at least if the bond/react can be used in my simulation. If not, I can continue using bond/break in my work.

Best Regards
Xu Hao

7

I cannot comment on fix bond/react. i’m copying the author, he may be able to help, when he finds the time. keep in mind, that it is currently christmas in the US and thus people tend to be with their family.

as for avoiding lost bond/atoms with fix bond/break. you have to figure out what is exactly happening. i.e. try to find a condition where you can cause lost bond atoms reliably and quickly after a known number of steps and then write out a trajectory frequently and visualize the positions.

for the most part you will need to avoid a very large change in force. this would usually be done by adjusting the bond breaking criterion. if you set the cutoff too large, then the change in force will be very large. another issue is the timestep. if the time step is too large, your system will not be able to handle the faster motion due to the bond breaking reaction correctly and you will get bad/diverging time integration which can lead to lost atoms/bond atoms. using a large communication cutoff and more frequent neighbor list updates can help, but in my personal opinion the suitable bond breaking criterion and a shorter timestep should be the core measures to cause stable simulations.

axel.

8

  1. Yes, bond/react is compatible with fix deform

  2. Correct, three reactions should be defined in that case. If C-C bond breaks (presumably a chain scission?), is it really possible for an adjacent O-C scission to occur? If so, it sounds like another reaction would have to be defined for that

9

Thanks so much Jacob and Axel, and sorry for disturbing you in holidays.

  1. Thanks for the information.

  2. You are right. In a chain scission only one bond in a -O-C-C-N- segment will break. As I see from the lammps example, the molecule templates include not only bonds but also angles, dihedrals, etc. So when I prepared the pre-reacted and reacted molecule templates for the O-C reaction, for example, will the templates also include the adjacent atom and bond (e.g., C-C) information ? If so, when a C-C bond is broken during tensile, will it influence the templates for the O-C ? I may misunderstood bond/react. If the templates for the three different reactions are totally independent from each other, my concerns will be not necessary.

Another question is that since there are thousands of atoms involved in the three reactions,should my templates include the topology information involving all the atoms ? If so there seems to be a tedious work. And I also see atom coordinates should be written in the templates. Since the atom coords continuously change during simulation, which state do the coords in the templates correspond to ?

Thanks a lot. If I can make sure that bond/react is suitable for my simulation, I can keep studying the details of this command anyway.

Best regards
Hao

10

It does sound like you are misunderstanding how bond/react works. Reaction templates only include atoms that are affected by the reaction (you may want to review the relevant parts of the doc page about this), and templates coordinates are not used by any current bond/react features (coords need not be included at all).

I am not sure i understand your first couple questions, but yes it sounds like these reactions are independent (and bond/react prevents overlapping reactions from occurring at the same time)

11

Many thanks for your explanation, Jacob. I will further study the details of this command.

Regards
Hao

12

Dear Jacob,

I think I am not more clear with the bond/react commmand. However, further concerns are that in the -O-C-C-N- chain that I introduced before, I understand three independent reactions should be specified to let one of the three bonds to break. However, there are about 1000 such chains in my system with equal possibility to break. So should I define 3*1000=3000 reactions, each including molecule templates and map files ? I guess so because in the map files the exact target IDs of the bonding atoms should be given, as I read from the lammps examples (i.e.,user/misc/bond_react). If so I am afraid the job amount will be too huge for me to use bond/react.

Regards
Hao

13

No, the reaction template atom IDs do not equal the simulation atom IDs. That would defeat the purpose of using templates.

Instead, just the atom types and connectivity of the pre-reaction template are used to find reaction sites in the simulation. Then, the simulation reaction site is updated to match the post-reaction template.

For example, the large nylon example user/misc/bond_react shows how a single polymerization reaction can build many long polymer chains.

Jake

14

Dear Jacob,

I have run some simulations using bond/react for the bond scissions in my system including the -O-C-C-N- chains. I have observed a sudden drop in the stress curve where should correspond to a bond breaking event (I am quite sure because by using bond/break command I got the same result, and one bond was deleted in the output bond information using compute/local command). However, the problem now is by using bond/react I didnt see any change of the output bond topology no matter by using compute/local or write_data. In other words the output information implies that no bond was broken. It is strange to me. Is there any reason for this ? For example, is bond/react using another way to calculate and output bond scission events ?

Regards
Hao

15

Hi Hao,

The outputted bond topology should also be updated, with bond/react. Are there a different number of bonds between your pre-reaction and post-reaction templates? If not, you may need to send a simple example for me to comment further.

Jake

16

Thanks Jacob, the bond numbers are different. And my input, template and map files are as follows:

units real
atom_style full
dimension 3
#newton on
boundary p p p

Style

pair_style hybrid buck/coul/long 15.0 lj/cut/coul/long 15.0 ##############

pair_modify mix arithmetic

bond_style hybrid harmonic morse
angle_style harmonic
dihedral_style charmm
improper_style umbrella
kspace_style pppm 1.0e-5
special_bonds dreiding

read_data relax_500ps_extra.lmps

pair_coeff * *…

neighbor 2.0 bin
neigh_modify delay 0 every 1 check no
timestep 0.5

variable tmp equal “lx”
variable L0 equal ${tmp}
variable strain equal “(lx-v_L0)/v_L0”
variable p1 equal “v_strain”
variable p2 equal “-pxx/101.01325"
variable p3 equal "-pyy/101.01325”
variable p4 equal “-pzz/10*1.01325”
variable p5 equal “lx”
variable p6 equal “ly”
variable p7 equal “lz”
variable p8 equal “temp”
variable t2 equal “epair”
variable t3 equal “ebond”
variable t4 equal “eangle”
variable t5 equal “edihed”
variable t6 equal “evdwl”
variable t7 equal “ecoul”
variable t8 equal “elong”

group bb type 2 3 8 9 # atom types involved in bond scission

molecule mol1 unreacted.data_template
molecule mol2 reacted.data_template

fix 1 all bond/react stabilization yes nvt_grp .03 &
react rxn1 all 2500 1.7625 4 mol1 mol2 map stabilize_steps 2000

fix 2 nvt_grp_REACT nvt temp 300.0 300.0 100.0

fix tctrl bond_react_MASTER_group temp/rescale 50 300 300 50 1

fix 3 all deform 1 x erate 0.000005 remap x units box #strain rate=5e9

thermo_style custom step temp etotal ke pe press lx ly lz
thermo 5000

compute 101 bb property/local btype batom1 batom2
dump 101 all local 10000 bond.dump index c_101[1] c_101[2] c_101[3]
fix def1 all print 700 “{p1} {p2} {p3} {p4} {p5} {p6} {p7} {p8}” file mech.txt screen no
fix def2 all print 700 “{p1} {t2} {t3} {t4} {t5} {t6} {t7} {t8}” file energy.txt screen no

run 1400000

write_data traj.data

unfix 1
unfix 2
unfix 3
uncompute 101
unfix def1
unfix def2

Molecule template before reaction:

this is a molecule template

26 atoms
25 bonds
43 angles
54 dihedrals
24 impropers

Types

1 2
2 9
3 3
4 2
5 2
6 5
7 6
8 6
9 6
10 6
11 6
12 7
13 2
14 2
15 5
16 2
17 6
18 6
19 6
20 1
21 1
22 1
23 1
24 1
25 6
26 6

Bonds

1 7 3 4
2 7 1 10
3 7 1 11
4 14 1 2
5 6 14 17
6 5 14 16
7 9 14 15
8 7 16 18
9 7 16 19
10 14 16 2
11 6 5 9
12 5 5 1
13 9 5 6
14 6 4 7
15 6 4 8
16 5 4 5
17 11 6 12
18 5 13 14
19 1 21 22
20 3 21 25
21 1 20 21
22 1 20 24
23 15 20 2
24 3 24 26
25 1 23 24

Angles

1 12 10 1 5
2 12 11 1 5
3 11 10 1 11
4 28 10 1 2
5 28 11 1 2
6 27 5 1 2
7 29 1 2 20
8 30 16 2 1
9 29 16 2 20
10 8 13 14 17
11 7 13 14 16
12 14 13 14 15
13 12 17 14 16
14 16 17 14 15
15 14 16 14 15
16 12 18 16 14
17 12 19 16 14
18 11 18 16 19
19 28 18 16 2
20 28 19 16 2
21 27 14 16 2
22 1 20 24 23
23 3 20 24 26
24 3 23 24 26
25 1 21 20 24
26 31 21 20 2
27 31 24 20 2
28 8 4 5 9
29 7 4 5 1
30 14 4 5 6
31 12 9 5 1
32 16 9 5 6
33 14 1 5 6
34 11 3 4 7
35 11 3 4 8
36 12 3 4 5
37 9 7 4 8
38 8 5 4 7
39 8 5 4 8
40 20 5 6 12
41 1 20 21 22
42 3 20 21 25
43 3 22 21 25

Dihedrals

1 20 10 1 5 6
2 20 11 1 5 6
3 44 10 1 2 20
4 44 11 1 2 20
5 42 5 1 2 20
6 43 5 1 2 16
7 45 10 1 2 16
8 45 11 1 2 16
9 12 13 14 16 18
10 12 13 14 16 19
11 18 17 14 16 18
12 18 17 14 16 19
13 39 13 14 16 2
14 40 17 14 16 2
15 41 15 14 16 2
16 20 18 16 14 15
17 20 19 16 14 15
18 45 18 16 2 1
19 44 18 16 2 20
20 45 19 16 2 1
21 44 19 16 2 20
22 43 14 16 2 1
23 42 14 16 2 20
24 47 23 24 20 2
25 48 26 24 20 2
26 1 24 20 21 22
27 2 24 20 21 25
28 1 21 20 24 23
29 2 21 20 24 26
30 46 21 20 2 1
31 46 24 20 2 1
32 46 21 20 2 16
33 46 24 20 2 16
34 12 1 5 4 7
35 12 1 5 4 8
36 12 4 5 1 10
37 12 4 5 1 11
38 18 9 5 1 10
39 18 9 5 1 11
40 24 4 5 6 12
41 25 9 5 6 12
42 24 1 5 6 12
43 39 4 5 1 2
44 40 9 5 1 2
45 41 6 5 1 2
46 15 3 4 5 9
47 30 3 4 5 1
48 17 3 4 5 6
49 18 7 4 5 9
50 20 7 4 5 6
51 18 8 4 5 9
52 20 8 4 5 6
53 47 22 21 20 2
54 48 25 21 20 2

Impropers

1 8 4 3 7 8
2 9 4 3 7 5
3 9 4 3 8 5
4 10 5 1 9 6
5 23 16 14 18 2
6 23 16 14 19 2
7 10 14 16 17 15
8 25 2 16 20 1
9 2 21 20 22 25
10 2 24 20 23 26
11 6 4 5 7 8
12 23 1 5 10 2
13 23 1 5 11 2
14 8 1 10 11 5
15 24 1 10 11 2
16 18 5 4 9 1
17 11 5 4 9 6
18 19 5 4 1 6
19 18 14 13 17 16
20 11 14 13 17 15
21 19 14 13 16 15
22 26 20 21 24 2
23 8 16 18 19 14
24 24 16 18 19 2

Molecule template after reaction:

this is a molecule template

26 atoms
24 bonds
38 angles
40 dihedrals
20 impropers

Types

1 2
2 9
3 3
4 2
5 2
6 5
7 6
8 6
9 6
10 6
11 6
12 7
13 2
14 2
15 5
16 2
17 6
18 6
19 6
20 1
21 1
22 1
23 1
24 1
25 6
26 6

Bonds

1 7 3 4
2 7 1 10
3 7 1 11
4 6 14 17
5 5 14 16
6 9 14 15
7 7 16 18
8 7 16 19
9 14 16 2
10 6 5 9
11 5 5 1
12 9 5 6
13 6 4 7
14 6 4 8
15 5 4 5
16 11 6 12
17 5 13 14
18 1 21 22
19 3 21 25
20 1 20 21
21 1 20 24
22 15 20 2
23 3 24 26
24 1 23 24

Angles

1 12 10 1 5
2 12 11 1 5
3 11 10 1 11
4 29 16 2 20
5 8 13 14 17
6 7 13 14 16
7 14 13 14 15
8 12 17 14 16
9 16 17 14 15
10 14 16 14 15
11 12 18 16 14
12 12 19 16 14
13 11 18 16 19
14 28 18 16 2
15 28 19 16 2
16 27 14 16 2
17 1 20 24 23
18 3 20 24 26
19 3 23 24 26
20 1 21 20 24
21 31 21 20 2
22 31 24 20 2
23 8 4 5 9
24 7 4 5 1
25 14 4 5 6
26 12 9 5 1
27 16 9 5 6
28 14 1 5 6
29 11 3 4 7
30 11 3 4 8
31 12 3 4 5
32 9 7 4 8
33 8 5 4 7
34 8 5 4 8
35 20 5 6 12
36 1 20 21 22
37 3 20 21 25
38 3 22 21 25

Dihedrals

1 20 10 1 5 6
2 20 11 1 5 6
3 12 13 14 16 18
4 12 13 14 16 19
5 18 17 14 16 18
6 18 17 14 16 19
7 39 13 14 16 2
8 40 17 14 16 2
9 41 15 14 16 2
10 20 18 16 14 15
11 20 19 16 14 15
12 44 18 16 2 20
13 44 19 16 2 20
14 42 14 16 2 20
15 47 23 24 20 2
16 48 26 24 20 2
17 1 24 20 21 22
18 2 24 20 21 25
19 1 21 20 24 23
20 2 21 20 24 26
21 46 21 20 2 16
22 46 24 20 2 16
23 12 1 5 4 7
24 12 1 5 4 8
25 12 4 5 1 10
26 12 4 5 1 11
27 18 9 5 1 10
28 18 9 5 1 11
29 24 4 5 6 12
30 25 9 5 6 12
31 24 1 5 6 12
32 15 3 4 5 9
33 30 3 4 5 1
34 17 3 4 5 6
35 18 7 4 5 9
36 20 7 4 5 6
37 18 8 4 5 9
38 20 8 4 5 6
39 47 22 21 20 2
40 48 25 21 20 2

Impropers

1 8 4 3 7 8
2 9 4 3 7 5
3 9 4 3 8 5
4 10 5 1 9 6
5 23 16 14 18 2
6 23 16 14 19 2
7 10 14 16 17 15
8 2 21 20 22 25
9 2 24 20 23 26
10 6 4 5 7 8
11 8 1 10 11 5
12 18 5 4 9 1
13 11 5 4 9 6
14 19 5 4 1 6
15 18 14 13 17 16
16 11 14 13 17 15
17 19 14 13 16 15
18 26 20 21 24 2
19 8 16 18 19 14
20 24 16 18 19 2

Map file:

this is a nominal superimpose file

4 edgeIDs
26 equivalences

BondingIDs

1
2

EdgeIDs

3
13
22
23

Equivalences

1 1
2 2
3 3
4 4
5 5
6 6
7 7
8 8
9 9
10 10
11 11
12 12
13 13
14 14
15 15
16 16
17 17
18 18
19 19
20 20
21 21
22 22
23 23
24 24
25 25
26 26

Thanks.

Hao

17

I won’t be able to test this simulation without the data file and full input file.