How to use read_restart along with loop

LAMMPS LAMMPS Beginners
, ,
1

Dear All,

I am providing an explanation of the problem, along with my program’s structure and sharing my input file for reference. I would greatly appreciate your assistance in incorporating the read_restart command and defining a loop that initiates from the equilibrium data file loaded by read_restart

Here is problem explanation:

I conducted a simulation of a nanotube using DPDe methodology. The interaction potential between bead-bead pairs within the nanotube’s wall can be adjusted by varying the bond coefficient. This coefficient serves as a pivotal parameter for manipulating the Young’s Modulus of the nanotube. To determine the Young’s Modulus, I subjected both ends of the CNT to opposing forces, initiating elongation.

However, to systematically explore the impact of varying bond coefficients on the nanotube’s mechanical properties, I need to loop over these coefficients. This loop, designated by the variable “loopa,” should encompass different bond coefficient values. Additionally, I must commence the elongation process from an equilibrium state, necessitating the utilization of the read_restart command.

My challenge lies in sequencing these steps effectively. I must ensure that the loop begins after all regions and computes are defined to prevent the “compute/region ID is reused” error. However, I also need the read_restart command at the outset to establish the initial state of the simulation. This predicament has left me uncertain about how to integrate the loop and read_restart command seamlessly to fulfill my objectives.

My program proceeds in the following steps:

  1. Initialization and Equilibration: Initially, the program reads the data, including bond coefficients, and conducts an equilibration phase for the carbon nanotube (CNT) at room temperature. This step ensures that the system reaches a stable equilibrium state before further manipulation.
  2. Elongation of the CNT: Subsequently, the CNT is subjected to elongation. This involves applying opposing forces to both ends of the nanotube to induce deformation. The response of the CNT to this elongation is monitored to understand its mechanical behavior under external loading conditions.
  3. Looping over Bond Coefficients: The program then enters a loop where it iterates over different bond coefficients. For each iteration, the CNT is elongated from its equilibrium position using the specific bond coefficient value corresponding to that iteration. This allows for a systematic exploration of how variations in bond coefficients influence the mechanical properties of the nanotube.

Inputfile:

variable T equal 1


 units 		        lj
 dimension	        3
 boundary	        p p p
 neighbor	        0.2 bin
 neigh_modify            every 1 delay 0 check yes
 neigh_modify            one 1500


 atom_style	       hybrid edpd full
 bond_style             harmonic
 angle_style            harmonic
 pair_style             edpd 1.58 9872598 



 # Reset box and simulation parameters

clear
region mybox  block -15 15 -15 15 -15 15 units box

read_restart restart.equi
Include parmcnt.lammps

               
 #GROUPING

group CNT    type 1


comm_modify vel yes cutoff 5


#RECENTER THE CNT

variable CNT_xcm equal -1*xcm(CNT,x)                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                       
variable CNT_ycm equal -1*xcm(CNT,y)
variable CNT_zcm equal -1*xcm(CNT,z)
displace_atoms CNT move ${CNT_xcm} ${CNT_ycm} ${CNT_zcm}

#CHANGE BOX BOUNDARIES

 change_box all x final -15 15 y final -15 15 z final -15 15

variable zmax equal bound(CNT,zmax)-1.5
variable zmin equal bound(CNT,zmin)+1.5
variable zmidmin equal ((${zmax}+${zmin})/2)-0.6
variable zmidmax equal ((${zmax}+${zmin})/2)+0.6
region rtop block INF INF INF INF ${zmax} INF
region rbot block INF INF INF INF INF ${zmin}
region rmid block INF INF INF INF ${zmin} ${zmax}
region zmid block INF INF INF INF ${zmidmin} ${zmidmax}


  group CNT_top region rtop
 group CNT_bot region rbot
 group CNT_mid region rmid
 group CNT_zmid region zmid


 set region rtop type 2
 set region rbot type 3
 set region zmid type 1


 #THE MOLECULAR DYNAMICS
velocity CNT_mid create ${T} 48455 mom yes rot yes
fix mynve all nve
compute Tmid CNT_mid temp

fix myber CNT_mid temp/berendsen ${T} ${T} 100
fix_modify myber temp Tmid


label       loopa
 variable    a loop 15
 variable    c equal ${a}
 variable    b equal 100*${a}
 bond_coeff  1 ${b} 0.77


  #RESTRAIN THE MOTION OF ATOMS AT THE EDGES
  #(FOR TWIST WE WANT TO FIX ONE EDGE ONLY AND TWIST THE OTHER        EDGE)

 fix mysf2 CNT_top setforce 0 0 0
 fix mysf3 CNT_bot setforce 0 0 0
 velocity CNT_top set 0 0 0
 velocity CNT_bot set 0 0 0


 variable L equal xcm(CNT_top,z)-xcm(CNT_bot,z)


 fix at2 all ave/time 10 100 1000 v_c v_j v_b v_f v_L file len${a}.txt
fix at3 all ave/time 10 100 1000 v_c v_j v_b v_f f_mysf2[3] f_mysf3[3]  file force${a}.txt


 variable mc equal 5*1.99*(10)^(-26)
  variable rcutoff equal 6.4633*10^(-10)
 variable kB equal 1.38*10^(-23) 
 variable tem equal 300
 variable epsilon equal (v_kB*v_tem)
 variable tau0 equal sqrt(v_mc*(v_rcutoff)^2/(v_epsilon))
 variable funit equal (v_mc*v_rcutoff)/(v_tau0)^2 
 variable cntA equal 1.23654*(10)^(-18)
 variable lini equal 22.1782*(v_rcutoff)

 variable stress equal ((abs(f_mysf2[3])+abs(f_mysf3[3]))*v_funit)/(v_cntA)
 variable strain equal (v_L*v_rcutoff-v_lini)/(v_lini)

 fix at4 all ave/time 1 1 1 v_b v_f v_stress  file stress${a}.txt
 fix at5 all ave/time 1 1 1 v_b v_f v_strain  file strain${a}.txt


 #CHECK THE TEMP OF NON-FROZZEN GROUP OVER TIME

  fix at1 all ave/time 10 100 1000 c_Tmid file output_temperature_middle_group.dat


  # RUN A SIMPLE EQULIBRATION TO BRING THE SYS TO THE REQUIRED TEMP    BEFORE APPLYING ANY DEFORMATION

  thermo 100
  thermo_modify temp Tmid

 timestep 0.01

 run 1000

 #Elongation


velocity CNT_top set 0 0 0.005
velocity CNT_bot set 0 0 -0.005


 run 100


 next        a
 clear
 jump        SELF loopa

I could not find any example in GitHub that people used loops and read_restart with same aim. If you know of any example that can help could please provide me with the link.

2

What has this got to do with it?

Your quoted input file is not very readable (didn’t you check it??) and that is because you didn’t quote it correct. The forum guidelines explain how to do it properly.

3

I used it to rewrord the question
ignore it please. I was editing it when you answered :slight_smile:

4

Hi Delaram,
If you can’t get it to work with LAMMPS internal loop system, there is always the possibility to launch the simulation with bash. Bash allows you to create for loops, and also to provide the values of lammps variables. This way, there wont be any issue of the type “compute/region ID is reused”.
Best,
Simon

1 Like
5

This can be dealt with rather easily. There is a variable function is_defined().
You can just define an equal style variable or use an immediate variable expansion to with is_defined(compute, someid) and it will return either 1 or 0 depending on whether that specific compute exists or not. The same applies to other categories like fixes, regions, dumps, groups, and variables.
Please note that with groups you may always call group someid delete without causing an error, even if the group “someid” does not exist.

Also for debugging, you can use the “info” command to see what is defined and what settings are made.

This is just standard computing programming stuff: You create a small test system that can run quickly and see what works empirically in terms of syntax for a non-functional small system with short trajectories. As mentioned above, you can use the info command to keep track of things and the variable function then to test for certain conditions in your script.

1 Like
6

Thanks a lot

7

Dr. akohlmey, The other error I face is : Bond_coeff command before simulation box is defined (src/input.cpp:1376)

The problem is I have to have the loop to start from the equilibrium state (restart.equi) on each step. that is the main reason I have used read_restart. I need the elongation to start from the equilibrium position, not the elongated position from the previous step. I think to do that, I have to place the loop just before read_restart restart.equi, but that doesn’t work. since the simulation box is defined in restart.equi. Even when I place the loop after the read_restart command the same error pops up. I was expecting the next step of the loop to jump to the first line of the loop, but it seems it starts from the beginning of the input file and that is why it produces the same error. Can you help me understand this please, and have the loop start from the equilibrium on each step?
I read examples and I read the documentation, asked people, but having a hard time make this work.

8

You didn’t read the documentation well enough. Please carefully re-read the documentation of the read_restart command. What you are apparently doing makes no sense considering what read_restart does. The error message is giving you the important hint of what you have to do.

Who said learning how to do simulations was easy?

9

Due to documentation I need to respecify the bond_coeffs in the restart file, but the my bond coeffs are changing due to the loop. I am not sure how to do that. Could you please give me a hint?

" Note that some force field styles (pair, bond, angle, etc) do not store their coefficient info in restart files. Typically these are many-body or tabulated potentials which read their parameters from separate files. In these cases you will need to re-specify thepair_coeff, bond_coeff, etc commands in your restart input script. The doc pages for individual force field styles mention if this is the case. This is also true of pair_style hybrid (bond hybrid, angle hybrid, etc) commands; they do not store coefficient info."

10

Dear Simon,

I just noticed that the bond_coeffs need to be re-written in the restart file . I found this answer to the similar question from Axel:

You could try one of the following:

  1. you could try using write_data and read_data instead of restart files.
  2. you could set up a minimal input script and then write out all XXX_coeff statements to a file with “write_coeff” and then using the “include” command to include those coeff settings after continuing from a restart file
  3. you can write a little python/perl/awk/bash/whatever script/program that processes the XXXCoeff sections of your data file into a file with XXX_coeff statements like in option 2). Most of them implicitly call the corresponding XXX_coeff command with the line of the data file appended. For the PairCoeff section there is only one atom type given, so the atom type has to be duplicated to match the requirements of the pair_coeff command and the mixed terms are automatically generated - like when using the data file - through the selected mixing rule. If you are using a PairIJCoeff section, that this section lists two atom types and has all mixed terms explicitly included, so just prepending pair_coeff would be sufficient.

Options 1) and 2) require that the used styles do support writing to a data file (correctly). This may depend on the LAMMPS version in use and the specific pair/bond/angle/dihedral/improper style.

since I am looping over bond_coeffs I am not sure how to do step 2. Writing a bash as you said is another option, which I have no idea how to start with. Is there any example or link to a tutorial or any further explanation of how it is done to learn from?

11

You are not thinking logically and you are not reading the documentation correctly.
You cannot set the bond coeffs before read_restart, so …?

What the read_restart docs say is that the restart file may contain the bond coeffs and will automatically re-apply them, but in which way does that interfere with you needing to change the bond coeffs?

12

  1. I have to loop over bond_coefficients in my simulation. 2. I need each loop to start from the equilibrium. 3. I can do that by creating a restart file in the equilibrium. 4. restart file contains the box information so I have to use it in the very beginning of the input file. 5. error is “bond_coeff command before simulation box is used”. 6. I guess this happens because in the next-step of the loop It reads the bond_coeff and then the restart file, that is why it’s producing that error. 7 ( I thought maybe I should try one of the ways you mentioned to have the bond_coeffs in the restart file, but don’t know how to do it as the bond_coeffs are in the form of a loop) 8. Maybe I am wrong about step “6”, as the next step of the loop due to the command " jump SELF loopa" should start from the the first line of the loop which is defined after the restart file.

  2. I find it complicated because I want the loop to start from equilibrium (restart file), but I don’t know how to do that as the restart file containing simulation box should be defined in advance. ( I cannot use the read and read_restart together as well)

You are right that I do not have the problem of " Note that some force field styles (pair, bond, angle, etc) do not store their coefficient info in restart files. Typically these are many-body or tabulated potentials which read their parameters from separate files. In these cases you will need to re-specify thepair_coeff, bond_coeff, etc commands in your restart input script. The doc pages for individual force field styles mention if this is the case. This is also true of pair_style hybrid (bond hybrid, angle hybrid, etc) commands; they do not store coefficient info." but I thought respecifying bond_coeffs maybe the solution to this, instead of changing the bond_coeffs using the loop outside of the restart file.

This was the way I thought. Maybe completely wrong. That is why I am asking to get help from the people who have more experience and developed better programming logic.

13

As I mentioned before, you are making a Mt. Everest out of a mole hill. If you cannot use the bond_coeff command(s) before the read_restart, you just use it after. What is so difficult about figuring that out?
It is plain simple logical thinking.

In fact, using bond_coeff after read_restart is required since the restart will reset the bond_coeff values from the system that you generated the restart with. So even if it would be possible to specify them before reading the restart, it will have the desired effect because they would be overridden.

1 Like
14

Thank you Dr.akohlmey, But that is what I did. I did define the loop in which I am changing the bond_coeffs and also include the parameter file, all after the read-restart as required. The whole Mt.Everst I face is because I do define the bond_coeffs after the restart file and I get the error “bond_coeff command before simulation box is used”.

Can you give a hint please?

15

This is not possible. There is no box when you issue this command. Did you perhaps delete the box with “clear”.

At any rate, this has become extremely off-topic now and has turned into a simple issue of your writing a complex input script despite having insufficient skills to debug it and monitor what is going on.

I am not your adviser or tutor and this forum is not the place to teach these skills.

Thus I am closing this topic now.

16