[lammps-users] illegal shake command

[pinky_09]

Hi,

I am trying to do the shear simulation of a bio-polymer. But, it gives an error message:

ERROR: Illegal fix shake command (…/fix_shake.cpp:155)

LAMMPS (22 Aug 2018)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (…/comm.cpp:87)
using 1 OpenMP thread(s) per MPI task
#MD simulation cellulose under shear

units real
dimension 3
boundary s s s

atom_style full
neighbor 2.5 bin
neigh_modify delay 5

pair_style lj/charmm/coul/charmm 8.0 10.0
bond_style harmonic
angle_style harmonic
dihedral_style harmonic
improper_style harmonic
read_data data.relax
orthogonal box = (5.65259 -58.6119 -45.9213) to (43.1264 21.408 46.6191)
2 by 6 by 8 MPI processor grid
WARNING: Pair style in data file differs from currently defined pair style (…/read_data.cpp:580)
reading atoms …
9066 atoms
reading velocities …
9066 velocities
scanning bonds …
4 = max bonds/atom
scanning angles …
6 = max angles/atom
scanning dihedrals …
17 = max dihedrals/atom
reading bonds …
9441 bonds
reading angles …
17398 angles
reading dihedrals …
32314 dihedrals
4 = max # of 1-2 neighbors
9 = max # of 1-3 neighbors
18 = max # of 1-4 neighbors
22 = max # of special neighbors

region lower block INF INF INF -31 INF INF
region upper block INF INF -8 INF INF INF
group lower region lower
3737 atoms in group lower
group upper region upper
3992 atoms in group upper
group boundary union lower upper
7729 atoms in group boundary
group mobile subtract all boundary
1337 atoms in group mobile

#set group lower type 6
#set group upper type 7

Assign original velocities to atoms

compute new all temp
velocity all create 298.1 487639 temp new

Set up ensemble

fix 1 all nvt temp 298.1 298.1 100.0
fix 2 all temp/rescale 10 298.1 298.1 0.01 1.0
unfix 2
fix 2 all shake 0.0001 10 100 21 4 6 8 10 12 14 18 47 31
ERROR: Illegal fix shake command (…/fix_shake.cpp:155)
Last command: fix 2 all shake 0.0001 10 100 21 4 6 8 10 12 14 18 47 31

Any help would be appreciated.

Thank you.

Sincerely,
Pinky

[Steve_Plimpton1]

If you look at the fix shake doc page and the examples at the top of it, they all
use the b/a/t/m keywords. Your command does not, hence it is incorrect.

Steve

[pinky_09]

Thanks for the tip!

[JWG1]

Dear all,
I try to simulate a rapid solidification process (cooling rate is about 1K/ps) in a melt pool with a thermal gradient of 10K/nm. The radius of the melt pool is about 40 nm. timestep 2 fs is used in the simulations, so the Tdamp is 100*0.002=0.2. Variable Theat is a function of the coordinates of each atom and cooling time. The corresponding code is as follows,
compute ke1 all ke/atom
fix 1 all nve
fix 2 meltpool langevin v_Theat 3500 0.2 123456
run 500000
unfix 2
unfix 1
I dumped the ke1, and found that the kinetic energy ( it is randomly distributed in the melt pool, which is messy after all), is not distributed according to my setting similar to Theat, so there is no directional solidification and columnar grains in the melt pool.
Any suggestions are appreciated very much.
Best regards,
Wugui Jiang

[akohlmey]

if you believe that the atom variable is not working correctly, please construct a minimal (small and fast to run) testcase that demonstrates it convincingly and post it here.
otherwise, there is likely something wrong with the parts of the input you are not showing, but it is not possible to give advice on something you cannot see or easily reproduce.

axel.

[JWG1]

Thank Axel,

My input code for debug is followed as,

units metal
boundary p p p
atom_style atomic

lattice fcc 4.04527
region box block -25 25 -25 25 0 2 units lattice
create_box 1 box
create_atoms 1 box
mass 1 26.981540000

variable tmpx0 equal (xlo+xhi)/2
variable x0 equal {tmpx0} variable tmpy0 equal yhi variable y0 equal {tmpy0}
variable tmpz0 equal (zlo+zhi)/2
variable z0 equal ${tmpz0}

region soregion block -25 25 -5 25 0 2 units lattice ###solidification region

group soregion region soregion
group other subtract all soregion

pair_style eam/alloy
pair_coeff * * Al1.eam.alloy Al

compute ke1 all ke/atom ###
variable temp1 atom c_ke1/(1.5*8.617343e-5)
dump 1 all custom 5000 00.dump.lammpstrj id type x y z c_ke1 v_temp1 #

compute tempmelt soregion temp

print “Begin Minimize******************************************************************************”

minimize 1.0e-25 1.0e-25 10000 10000 #
undump 1
reset_timestep 0

timestep 0.002

velocity all create 293.0 451547 rot yes dist gaussian
thermo_style custom step temp pe ke etotal vol lx ly lz press pxx pyy pzz v_x0 v_y0 v_z0 v_tmpx0 v_tmpy0 v_tmpz0
thermo_modify lost ignore flush yes
thermo 1000

dump 2 all custom 5000 01.dump.lammpstrj id type x y z c_ke1 v_temp1 ###

print “Begin NVT 293 K ******************************************************************************”
fix 3 all nvt temp 293 293 0.2 drag 0.5 #
run 10000
unfix 3
undump 2
reset_timestep 0

change_box all boundary p s p
print “Begin heating ******************************************************************************”
fix h1 all nve
variable tmp equal “time” ###
variable tc equal {tmp} ### variable tt equal "time-v_tc" ###tt variable RR atom -y+v_y0 ### distance from the top variable Tr1 atom 2000-5*v_RR ### thermal gradient 5 K/A variable Tr atom (v_RR<=160)*v_Tr1+(v_RR>160)*293 ### 160 A variable Theat atom (v_tt>=100)*v_Tr+(v_tt<100)*(293+(v_RR<=160)*(v_Tr-293)*v_tt/100) ### 100 ps finish the heat process, 20 ps hold on thermo_style custom step temp c_tempmelt pe ke etotal vol lx ly lz press pxx pyy pzz density enthalpy v_tt v_x0 v_y0 v_z0 v_tmpx0 v_tmpy0 v_tmpz0 thermo_modify lost ignore flush yes thermo 1000 dump 3a all custom 5000 02.dump.lammpstrj id type x y z c_ke1 v_RR v_Tr v_Theat v_temp1 fix 4 other langevin 293 293 0.2 123456 fix 5 soregion langevin v_Theat 293 0.2 123456 run 60000 # {Nheat}
unfix 4
unfix 5
unfix h1
undump 3a

reset_timestep 0
print "Begin Cooling ***************************************************************************"
fix c1 all nve
variable tmp1 equal “time” ###
variable tcc equal ${tmp1} ###
variable tt1 equal “time-v_tcc” ###
variable Tcool atom ((v_Tr-v_tt11)<293)293+((v_Tr-v_tt11)>=293)(v_Tr-v_tt11) # 1K/ps cooling rate

thermo_style custom step temp c_tempmelt pe ke etotal vol lx ly lz press pxx pyy pzz density enthalpy v_tt1 v_x0 v_y0 v_z0 v_tmpx0 v_tmpy0 v_tmpz0
thermo_modify lost ignore flush yes
thermo 1000
dump 4a all custom 10000 03.dump.lammpstrj id type x y z c_ke1 v_Tcool v_Tr v_RR v_temp1

fix 5 other langevin 293 293 0.2 123456
fix 6 soregion langevin v_Tcool 293 0.2 123456
run 900000 # ${Ncool}

unfix 6
unfix 5

unfix c1
undump 4a

在 2021-01-20 12:20:33，“Axel Kohlmeyer” [email protected] 写道：

[JWG1]

Dear all,
any suggestions for this code?
thanks very much.

wugui jiang

发自我的华为手机

发件人： JWG [email protected]
日期： 2021年1月20日周三 晚上8:22
收件人： Axel Kohlmeyer [email protected]
抄送： LAMMPS Users Mailing List [email protected]
主 题： Re: [lammps-users] Kinetic energy of atom during fix langevin

[akohlmey]

It is not likely that somebody is going to debug your input for you.

axel.

[akohlmey]

this is not what I have asked for. I had asked for a convincing demonstration with a specific small and fast to run input that fix langevin with an atom style variable does not work the way it is supposed to work.
as a matter of demonstration please see the slightly modified input of the melt example and the two attached images.

this input assigns a target temperature with an atom style variable following a sine curve with two periods in y-direction.
the first image (image.0.jpg) shows the system in the initial state and the second image (image.2500.jpg) at the end of the (very short) simulation.
the color of the atoms represents the kinetic energy of each atom. you can clearly see the random (initial) distribution at step 0 and how this evolves over the simulation to two zones of high temperature (and disordered structure) as well as two zones of low temperature and regular structure. this corresponds to the prescribed input and cannot be explained other then by the langevin fix working exactly as expected.

the logical conclusion is that either how you compute your target temperature is not correct, or your expectation of how that target temperature would be represented in the simulation is not accurate. in either case, this is something for you to figure out and debug. hopefully, I have given you an example of how you can approach this.

good luck,
axel.

here is the input:

region box block 0 10 0 30 0 10
create_box 1 box
create_atoms 1 box
mass 1 1.0

velocity all create 3.0 87287

pair_style lj/cut 2.5
pair_coeff 1 1 1.0 1.0

neighbor 0.3 bin
neigh_modify every 10 delay 0 check no

fix 1 all nve

variable temp atom (sin(4PIy/ly)+1.0)
fix 2 all langevin v_temp 1.0 0.1 62342614
compute 1 all ke/atom
variable tatm atom (c_1/3.0)

dump 1 all image 2500 image.*.jpg v_tatm type view 60 -30 zoom 1.2 size 1024 768 ssao yes 623425 0.3

thermo 500
run 2500

[akohlmey]

[…]

I think my code is correct, so the result is also correct.

sorry, but this way of thinking makes it almost certain that you are making mistakes and are producing results that you don’t understand and cannot properly explain.

Maybe the result is sensitive to the system size and timestep.

speculation is no way to address such a problem. you should simplify your system (there are a lot of details in your input needlessly complex and convoluted and difficult to read) start with a simpler system and approach to your cooling experiment and then gradually make it more complex while checking at each step whether both you input and the effect is like you expect it. only then can you be certain that you have a correct input and you will have understood why you are not getting the desired results in case that does not happen.

in fact, I would claim that the more experienced somebody is as a researcher, the more paranoid that person will become about making mistakes. as nobody makes no mistakes and often you wonder why you didn’t notice even the most trivial oversights.

axel.

[JWG1]

Dear Axel,
thanks for your kindness. I think my code is correct, so the result is also correct. Maybe the result is sensitive to the system size and timestep.

Best,
wugui jiang

在 2021-01-21 09:17:08，“Axel Kohlmeyer” [email protected] 写道：

[JWG1]

Dear Axel,

Thank you for your reminders and suggestions. We will examine the script line by line. Thanks again.

在 2021-01-21 11:13:29，“Axel Kohlmeyer” [email protected] 写道：