coreshell minimisation error

LAMMPS LAMMPS Mailing List Mirror
1

Dear Users and Developers,

I am trying to run a coreshell system, with the Li species within the system as rigid ion. I have made sure the correct modules are loaded, and tested the example coreshell model provided, which runs fine. My data file is in the correct format, where the element, atom, and bond groups/types are defined as:

Element type 1=Li, 2=Ni, 3=Ni_shell, 4=Co, 5=Co_shell, 6=Mn, 7=Mn_shell, 8=O, 9=O_shell

Group type 1=Li, 2= Ni and Ni_shell, 3= Co and Co_shell, 4= Mn and Mn_shell, 5= O and O_shell

Bond type 1 = Ni and Ni_shell, 2=Co and Co_shell, 3= Mn and Mn_shell, 4= O and O_shell​

Where the bond between each individual core-shell has an individual bond number, such as:

Bonds

1 1 961 962
2 1 963 964
3 1 965 966

4 1 967 968
… … … …

768 1 2495 2496
769 2 2497 2498
770 2 2499 2500 etc.

I am fairly confident that this input is correct. The system does run, however doesn’t get past step 0 of the minimisation. Below is the output given, where the program doesn’t stop running or exit, but will not go pass this point. I have to manually stop the process.

Setting up cg style minimization …
Unit style : metal
Current step : 0
Per MPI rank memory allocation (min/avg/max) = 69.21 | 69.21 | 69.21 Mbytes
Step TotEng Temp Enthalpy Fmax Lx Ly Lz Volume Press
0 -nan 0 -nan nan 45.623566 24.789324 27.811644 31454.339 -nan

I believe the issue may lie in the parameter input file, however I can not see where the error is, nor can my colleagues. I have tried varying everything I can think off but I cannot find anything which doesn’t result in the above. I have this exact system and potentials working in dlpoly4 without any issue, so the potentials are not the cause. Any assistance would be very much appriciated.

#------------Variables and cell--------------------------#
clear
units metal #eV,atomic charge,angstroms,ps,kelvin,bars,g/mol
dimension 3
boundary p p p
atom_style full
fix csinfo all property/atom i_CSID
read_data coordinates.lammps fix NULL CS-Info #File name read in containing coordinates etc.
group cores type 1 2 4 6 8
group shells type 3 5 7 9
group lithium type 1
neighbor 2.0 bin
comm_modify vel yes
variable T1 equal 1 #Sets temperature of the system
variable Timer equal step*dt

#---------- Pair styles and electrostatics--------------#
pair_style buck/coul/long/cs 8.0 #A, rho, C
pair_coeff * * 0.0 1.000 0.00

pair_coeff 1 9 632.1018 0.2906 0.0000

pair_coeff 3 9 1582.5000 0.2882 0.00000
pair_coeff 5 9 1397.6300 0.3211 0.00000
pair_coeff 7 9 1329.8200 0.3087 0.00000
pair_coeff 9 9 22764.0000 0.1490 65.0000

bond_style harmonic
bond_coeff 1 93.7 0.0
bond_coeff 2 196.3 0.0
bond_coeff 3 95.0 0.0
bond_coeff 4 65.0 0.0

kspace_style ewald 1e-03 #Ewald precision

-------------- Run Minimization ---------------------#

reset_timestep 0
timestep 0.5
thermo 10 #prints data every x steps
thermo_style custom step enthalpy fmax lx ly lz vol press
min_style cg
minimize 1e-25 1e-25 5000 10000

​Many thanks,

Lucy

2

Dear Lucy,

NaN as enthalpy or force strongly suggests bad initial geometry, i.e. at least two atoms (cores/shells not bonded to one another) with exactly the same coordinates.
I suggest dumping atom energies and checking whether all of them are NaNs or just some. If the former is true, it is an indication of a systematic error, for example wrong IDs in the bond section - in LAMMPS they start with 1, I’m not sure what convention is used in dlpoly4.
The latter case means the problem is local, for example due to periodic boundary conditions.
I can’t really give any more specific advice without a complete input and data file.

One other thing - you haven’t stated which version of LAMMPS you use. If it isn’t the newest one, try updating it.

Michał

3

Dear Michał​,

Thank you for your response. I thought the same thing and have gone through each individual pair for the whole system, which seems fine and matches up to the pairs used in the converted dlpoly4 system. I have tried also to dump the individual energy values of the atoms though I have not done this before and the specific command does not seem to work yet.

The bonds IDs do start at 1, so I can not see a specific issue here. Dlpoly4 uses the atom ID rather than a number to identify.

I do currently have the latest version of LAMMPS, and have also tried older versions and on set ups where coreshell models have worked before (to check the correct packages are loaded).

I believe this is probably a silly minor mistake I have made somewhere, though I can’t find it, or there is some issue with mixing the rigid Li ion with the coreshell parameters for the rest of the system.

I have attached the input files, if anyone is able to take a look and see any glaring mistakes. I have made a lot of changes to the parameter.in file so this is simplified for error testing purposes.

Many thanks,

Lucy

parameters.in (2.96 KB)

coordinates.lammps (574 KB)

4

Dear Michał,

Thank you for your response. I thought the same thing and have gone through each individual pair for the whole system, which seems fine and matches up to the pairs used in the converted dlpoly4 system. I have tried also to dump the individual energy values of the atoms though I have not done this before and the specific command does not seem to work yet.

The bonds IDs do start at 1, so I can not see a specific issue here. Dlpoly4 uses the atom ID rather than a number to identify.

I do currently have the latest version of LAMMPS, and have also tried older versions and on set ups where coreshell models have worked before (to check the correct packages are loaded).

I believe this is probably a silly minor mistake I have made somewhere, though I can't find it, or there is some issue with mixing the rigid Li ion with the coreshell parameters for the rest of the system.

I have attached the input files, if anyone is able to take a look and see any glaring mistakes. I have made a lot of changes to the parameter.in file so this is simplified for error testing purposes.

lucy,

thanks for providing a well laid out, readable, and complete input
deck. this helps so much in checking things out. i wish more people
here would follow your example.
the most obvious issue is, that you are using the wrong pair style.
the regular coulomb pair styles cannot handle coreshell systems. thus
you need to replace buck/coul/long with buck/coul/long/cs

also, your ewald convergence is too loose (why not pppm, which should
be much faster?), 1.0e-4 is barely acceptable for constant volume, but
for variable volume, you need better converged stresses and that
requires more like 1.0-6. and your time step seems very large (keep in
mind, that for metal units the unit of time is 1 ps)

axel.

5

Dear Axel,

thanks for providing a well laid out, readable, and complete input deck. this helps so much in checking things out. i wish more people here would follow your example.

Thank you, I understand how important it is to give as much information as possible to help resolve these sorts of issues.

the most obvious issue is, that you are using the wrong pair style. The regular coulomb pair styles cannot handle coreshell systems. thus you need to replace buck/coul/long with buck/coul/long/cs

This is excellent, thank you. I had look at this pair style originally though the potentials I have don't state a 'D' term, and therefore I chose the pair style which fit with the potential terms I have. I have looked through the lammps documentation and am trying to work out what this term is so I can calculate it. I understand there are two additional terms in this pair potential, the first of value 0.00 and the second, the cs term, which is additional to A, rho, and C. I can not seem to find anywhere in the documentation what this value D is, and the values in the example do not match any other data expressed. I may be being entirely clueless on this point, but would you be able to clarify what this term is for me please?

also, your ewald convergence is too loose (why not pppm, which should be much faster?), 1.0e-4 is barely acceptable for constant volume, but for variable volume, you need better converged stresses and that requires more like 1.0-6. and your time step seems very large (keep in mind, that for metal units the unit of time is 1 ps)

This is just one of the parameter files I used, varying all the terms I could trying to get the system to work. The original parameters did include pppm to 1^-6 and a time step closer to 0.001 ps (also 0.0001 ps). As I have used ewald in dlpoly previously to 1^-3 and 1^-6 I attempted to vary these also.

Many thanks,
Lucy