Using Hybrid MD + MC for polydisperse system with the help of fix/atom swap LAMMPS (29 Sep 2021)

Dear LAMMPS Users and Developers,

I am trying to perform Hybrid MD + MC for polydisperse system. For this I firstly include a parameter file in which I defined the pair_coeff for each particle. Then in a loop I ran MD simulations for 100 steps and then 100 steps for MC swap, where for the MC swap I used random numbers in the type command which lies in the range of total number of atoms. But if the two random numbers are same then swap will accept with probability 1. So to keep these random numbers different is essential. I am trying different to apply this but I am not able to figure this out.

Also firstly equilibrating the system, I wanted to do production where I want the output according to a particular file. Earlier when I used only MD this can be done easily but not since I am also using MC swap I am encountering problem doing this as well.

I would be grateful if anyone can suggest something about this.

The sample script is a as follows :-

# 3d Lennard-Jones melt

package intel 	0
units		lj
atom_style	atomic
dimension	3
boundary	p p p

lattice 	sc 1.0

region		cube block 0 5 0 5 0 5 
create_box	125 cube
create_atoms	1 box


mass		* 1.0

label 		loop
variable 	i loop 125
set 		atom ${i} type ${i}
next 		i 
jump 		in.twoatoms loop

velocity	all create 0.30 561564 mom yes rot yes

timestep	0.001

#global cutoff
pair_style      lj/cut 1.25
include   	parameter.dat

neighbor        1.60 bin 

neigh_modify    every 1 delay 0 check yes

variable 	random1 equal random(1,125,1564848)
variable 	random2 equal round(v_random1)
variable 	random3 equal random(1,125,514447)
variable 	random4 equal round(v_random3)

variable	t equal temp


label		loopf
variable	k loop 1

fix 		1 all nvt/intel temp 0.30 0.30 0.1
run		100
unfix		1

label		loopswap
variable	j loop 100

#if		"${random3} >= ${random4}" then "jump in.twoatoms reassign_random"
fix 		2 all atom/swap 1 1 29494 0.30 types ${random2} ${random4}
#if		"random2 == {random4}" then "fix 2 all atom/swap 1 1 29494 0.30 types ${random2} ${random4}"
run		1
unfix		2

#print		"${random2} ${random4}" append info.dat

next		j
jump		in.twoatoms loopswap
next		k

jump		in.twoatoms loopf


reset_timestep 	0

variable 	tmsd file sort.txt
variable 	tms equal next(tmsd)

dump		1 all custom 1 position.txt id type xu yu zu
dump_modify	1 sort id every v_tms #first yes


label		loopf2
variable	k loop 5

fix 		1 all nvt/intel temp 0.30 0.30 0.1

#variable	l equal "round((v_k - 1 ) * 100)" 
#reset_timestep		${l}

#print		"${tms} " file info.dat

run		100

unfix		1

label		loopswap2
variable	j loop 100

fix 		2 all atom/swap 1 1 29494 0.30 types ${random2} ${random4}
run		1
unfix		2

#print		"${random2} ${random4}" append info.dat

next		j
jump		in.twoatoms loopswap2

next		k
jump		in.twoatoms loopf2

undump		1

Here the “parameter.dat” file includes the parameters and “sort.txt” the timesteps I want my output at.

Best Regards
Puneet

1 Like

First off, to avoid having the forum software incorrectly format your inputs, you need to quote it into triple backquotes “```” before and after the text block (like a python string block).

There are multiple issues with your random number generation. When defining an equal style variable you effectively are creating something that more resembles a function than a variable. So if you call a function like random() it will generate a new random number every time you evaluate the variable, so it cannot work the way you have set it up. You need to change the way you loop and the way how you define the variables. The following block of input will generate a pair of different random numbers within the given range.

label next_rand
variable type_a equal $(round(random(1,125,12343123)))
variable type_b equal $(round(random(1,125,43234256)))
if "${type_a} == ${type_b}" then "jump SELF next_rand"

Note that by using a $(…) expression the random number is not created when calling the equal style variables but rather whenever the input line is processed.

As for the dump. You are redefining the variable and the dump instance in every loop block. That will start it all over again all the time. Thus you need to define those outside the loop and you must not use the reset_timestep command inside your loop, either.

1 Like

Dear Alex,

Thank you for your help. I will change how the random numbers are defined. As for the dump command, I have tried many different types and also the method as suggested by you but I have still not figured it out. Since the particles positions are incremented only in MD, so I want my positions at the MD steps, but I have to use “run” command in “fix atom/swap” as well, which creates a problem.

I would be grateful if you could show in a block or script format how I can overcome the problem of dumping the positions at MD timesteps that are in a file named “sort.txt”.

Thanks and Regards
Puneet

Using fix atom/swap without also having time integration makes no sense.

Dear Alex,

What I mean to say is that the particle positions are not changed after the application of swap move only there radius changes.

Best Regards
Puneet

That is what I understood. It still does not make sense. You have to orchestrate your swaps differently, either through direct swapping with script commands or by using a modified version of the fix atom/swap that is implementing your swap algorithm. Due to your very unusual use of atom types, the logic implemented in the existing fix is not suitable and your use of the fix is some kind of “abuse” and thus not applicable.

Dear Alex,

Thanks again for your help. I will look into this.

Best Regards
Puneet

Dear Alex,

Can you suggest some correct way I can implement the method that I am trying to do in LAMMPS.

Best Regards
Puneet

I already did.

Yeah.

Dear Alex,

I want to ask one more question.

The LAMMPS document says for the “run N” command that “A value of N = 0 is acceptable; only the thermodynamics of the system are computed and printed without taking a timestep”, so if I use

fix 2 all atom/swap 1 1 29494 300.0 types 1 2
run 0
unfix 2

will the swap attempt happen, since the atoms are not displaced at all by this command.

Best Regards
Puneet

This is a superfluous question. You can answer this question yourself by monitoring the coordinates.

But as I mentioned before, this kind of use is “abuse” of the fix command which is supposed to be used as part of the run and not hacked to be used like a standalone command and thus very problematic. Whether you do a run 1 or a run 0 does not make a difference in that regard.

Dear Alex,

First of all I want to say that what you are doing is really a great work. The time you put in for the research community to grow, is extraordinary.

Secondly thanks for the idea to check about the “run 0” command.

Thanks and Regards
Puneet

This is how research works! You start with a hypothesis, conduct a suitable experiment, observe carefully, and either confirm your hypothesis, reject it, or construct a more specific experiment to resolve inconclusive observations. It makes no difference whether this applies to physical experiments, computer simulations, or thought experiments.

1 Like

Yeah. :smiley: