Help about Initial Structure optimization and Selection (LixMn2O4 MD Simulation)

The unit cell of LiMn2O4 considered in my simulations contains 56 ions (8 lithium atoms, 16 manganese atom and 32 oxygen atoms)

In the LixMn2O4 system, the Li-ions are positioned on the 8a tetrahedral sites, and the manganese ions are located on the 16d octahedral sites. All possible initial configurations with a combination of Li+, Mn3+ and Mn4+ are constructed. For example,for x = 0.125, the number of Li+ is one and the number of available sites for Li+ is eight. For the manganese ions, one Mn3+ and fifteen Mn4+s occupy sixteen sites. Thus, the total number of crystal configurations is 128. In this way the number of cases for x = 0.250, 0.375, 0.500, 0.625, 0.750, 0.875 and 1.00 is 3360, 31360, 127400, 244608, 224224, 91520 and 12870, respectively.

Here a large number of cases for each x value are resulted from the intrinsic MD methodology that treats Li+, Mn3+ and Mn4+ as different atoms. Since it is almost impossible to check the equilibrium energy of all these cases, I want calculate the potential energy of these initial configurations by running one time step of MD simulations and sort the cases according to energy levels.

Main questions is :

How can I have an input script for automatically change initial structures using a loop or etc?

How can i produce all initial configurations and test them energetically?

# MD Run LiMn2O4 input script

#------------------------------- Initialization -----------------------------------------

Clear

units metal

dimension 3

atom_style charge

boundary p p p

#------------------------Atomic Structure Definition---------------------------------------------

variable alat equal 8.24

lattice custom ${alat} &

basis 0.0000 0.0000 0.0000 basis 0.2500 0.2500 0.2500 basis 0.0000 0.5000 0.5000 basis 0.5000 0.0000 0.5000 &

basis 0.5000 0.5000 0.0000 basis 0.2500 0.7500 0.7500 basis 0.7500 0.2500 0.7500 basis 0.7500 0.7500 0.2500 &

basis 0.6250 0.6250 0.6250 basis 0.6250 0.1250 0.1250 basis 0.1250 0.6250 0.1250 basis 0.1250 0.1250 0.6250 &

basis 0.6250 0.8750 0.8750 basis 0.6250 0.3750 0.3750 basis 0.1250 0.8750 0.3750 basis 0.1250 0.3750 0.8750 &

basis 0.8750 0.6250 0.8750 basis 0.8750 0.1250 0.3750 basis 0.3750 0.6250 0.3750 basis 0.3750 0.1250 0.8750 &

basis 0.8750 0.8750 0.6250 basis 0.8750 0.3750 0.1250 basis 0.3750 0.8750 0.1250 basis 0.3750 0.3750 0.6250 &

basis 0.3750 0.3750 0.3750 basis 0.3750 0.8750 0.8750 basis 0.8750 0.3750 0.8750 basis 0.8750 0.8750 0.3750 &

basis 0.3750 0.6250 0.6250 basis 0.3750 0.1250 0.1250 basis 0.8750 0.6250 0.1250 basis 0.8750 0.1250 0.6250 &

basis 0.6250 0.3750 0.6250 basis 0.6250 0.8750 0.1250 basis 0.1250 0.3750 0.1250 basis 0.1250 0.8750 0.6250 &

basis 0.6250 0.6250 0.3750 basis 0.6250 0.1250 0.8750 basis 0.1250 0.6250 0.8750 &basis 0.1250 0.1250 0.3750 &

basis 0.8750 0.8750 0.8750 basis 0.8750 0.3750 0.3750 basis 0.3750 0.8750 0.3750 basis 0.3750 0.3750 0.8750 &

basis 0.6250 0.6250 0.8750 basis 0.6250 0.1250 0.3750 basis 0.1250 0.6250 0.3750 basis 0.1250 0.1250 0.8750 &

basis 0.6250 0.8750 0.6250 basis 0.6250 0.3750 0.1250 basis 0.1250 0.8750 0.1250 basis 0.1250 0.3750 0.6250 &

basis 0.8750 0.6250 0.6250 basis 0.8750 0.1250 0.1250 basis 0.3750 0.6250 0.1250 basis 0.3750 0.1250 0.6250

region LMObox block 0 2 0 2 0 2

create_box 4 LMObox

create_atoms 4 box basis 1 1 basis 2 1 basis 3 1 basis 4 1 basis 5 1 basis 6 1 basis 7 1 basis 8 1 &

basis 9 2 basis 10 2 basis 11 2 basis 12 2 basis 13 2 basis 14 2 basis 15 2 basis 16 2 &

basis 17 3 basis 18 3 basis 19 3 basis 20 3 basis 21 3 basis 22 3 basis 23 3 basis 24 3 &

basis 25 4 basis 26 4 basis 27 4 basis 28 4 basis 29 4 basis 30 4 basis 31 4 basis 32 4 &

basis 33 4 basis 34 4 basis 35 4 basis 36 4 basis 37 4 basis 38 4 basis 39 4 basis 40 4 &

basis 41 4 basis 42 4 basis 43 4 basis 44 4 basis 45 4 basis 46 4 basis 47 4 basis 48 4 &

basis 49 4 basis 50 4 basis 51 4 basis 52 4 basis 53 4 basis 54 4 basis 55 4 basis 56 4

mass 1 6.94

mass 2 54.94

mass 3 54.94

mass 4 15.999

group Lithium type 1

group Manganese3 type 2

group Manganese4 type 3

group Oxygen type 4

set group Lithium charge +1

set group Manganese3 charge +1.4

set group Manganese4 charge +2.4

set group Oxygen charge -1.2

# -----------------------------------Force Field-----------------------------------------

pair_style born/coul/long 7.5 15

pair_modify tail yes

kspace_style ewald 1e-05

pair_coeff 1 1 0.2637055573 0.317 2.34 1.0485818488 -0.4993247392

pair_coeff 1 2 0.2637055573 0.317 2.34 1.0485818488 -0.4993247392

pair_coeff 1 3 0.2637055573 0.317 2.34 1.0485818488 -0.4993247392

pair_coeff 1 4 0.2637055573 0.317 2.34 1.0485818488 -0.4993247392

pair_coeff 2 2 0.2637055573 0.317 2.34 1.0485818488 -0.4993247392

pair_coeff 2 3 0.2637055573 0.317 2.34 1.0485818488 -0.4993247392

pair_coeff 2 4 0.2637055573 0.317 2.34 1.0485818488 -0.4993247392

pair_coeff 3 3 0.2637055573 0.317 2.34 1.0485818488 -0.4993247392

pair_coeff 3 4 0.2637055573 0.317 2.34 1.0485818488 -0.4993247392

pair_coeff 4 4 0.2637055573 0.317 2.34 1.0485818488 -0.4993247392

neighbor .3 bin

neigh_modify every 10 delay 10 check yes page 1000000

# -----------------------------------Initial velocities-------------------------------

compute mytemp all temp

compute mypress all pressure mytemp

velocity all create 300 6767888 temp mytemp

# ----------------------------------------Running-------------------------------------------

###### Equilibrate Phase #####################################

fix 1 all npt temp 300 300 .01 iso 0 0 .1

fix_modify 1 press mypress

fix 2 all temp/rescale 10 300.0 300.0 10.0 1.0

fix_modify 2 temp mytemp

#Set boundary conditions to be stress-free

fix 3 all box/relax iso 0.0

thermo 100

thermo_style custom step temp etotal pe ke press lx

thermo_modify temp mytemp press mypress lost error

timestep 0.001

dump 1 all atom 1 dump_minimize.lammpstrj

run 10000

It would be easier to write a shell script to supply the variables to a template input script, then use the same shell script to drive the LAMMPS simulations.

If you just want single point energy, why do you have a “run 10000” at the end?

Ray

It would be easier to write a shell script to supply the variables to a
template input script, then use the same shell script to drive the LAMMPS
simulations.

...or use python and the python interface to LAMMPS?

in general, one should always keep in mind that despite its
flexibility, the LAMMPS scripting language is not a general purpose
script language and thus has limited applicability to all such complex
scripting tasks and a proper scripting language is much more adequate.

axel.

Thanks for your reply

Because for all initial configurations(about 700,000 initial configurations), equilibrium energy is not accessible I want just single point energy and 1 MD step is correct. Now I have 3 questions :

What is the shell script?

Can I use MATLAB interface instead of Python interface for my problem?

And what is the difference between these two methods?

Thanks alot

I suggest you use any programming language you are comfortable
with and write a small program that generates the 700,000 LAMMPS input
files and runs them, one at a time, and renames the produced log.lammps
file to log.N, where N runs from 1 to 700,000.

Then you can post-process the log files however you wish, e.g. with another
small program.

Python is ideal for this. You don’t need the Python interface to LAMMPS,
just Python. But if you know Matlab better, I imagine you
can do it in Matlab as well.

In Python it would be as simple as something like this:

for i in xrange(700000):
infile = "in.d" i
fp = open(infile,“w”)

… code to write an input script into infile

cmd = "lmp_g++ < s" infile
out = commands.getoutput(cmd)
os.rename(“log.lammps”,infile)

Steve