Potential energy equal zero

LAMMPS LAMMPS Mailing List Mirror
1

Hi,

I am working towards developing interatomic potentials for certain alloys and writing my par_styles in c++. The code that I wrote is a mix of c++ and python. For model validation, I have written energy and force computations in python, which are returned to pair_style written in c++.

The issue is that pair_style works fine for simple calculations such as cohesive energy and lattice constant. The correct energy is returned from python code as shown below

('ENERGY IN PYTHON ', array([[-1235.68522697]]))

Per MPI rank memory allocation (min/avg/max) = 4.004 | 4.004 | 4.004 Mbytes

Step Temp PotEng KinEng Lx Ly Lz Press Pxx Pyy Pzz c_eatoms

0 0 -1235.6852 0 49.2 8.5216899 24.6 0 0 0 0 -1235.6852

But it becomes zero, when I use the same potential for something like phonon dispersions, as shown below

('ENERGY IN PYTHON ', array([[-1409.1978995]]))

Per MPI rank memory allocation (min/avg/max) = 3.919 | 3.919 | 3.919 Mbytes

Step Temp PotEng KinEng Press Pxx Pyy

0 300 0 7.7168307 748.18204 776.72401 758.39058

I am not sure why is it happening. Can you please help me in resolving this issue?

Regards,

Gurjot Dhaliwal

2

Hi,

I am working towards developing interatomic potentials for certain alloys and writing my par_styles in c++. The code that I wrote is a mix of c++ and python. For model validation, I have written energy and force computations in python, which are returned to pair_style written in c++.

The issue is that pair_style works fine for simple calculations such as cohesive energy and lattice constant. The correct energy is returned from python code as shown below

('ENERGY IN PYTHON ', array([[-1235.68522697]]))

Per MPI rank memory allocation (min/avg/max) = 4.004 | 4.004 | 4.004 Mbytes

Step Temp PotEng KinEng Lx Ly Lz Press Pxx Pyy Pzz c_eatoms

0 0 -1235.6852 0 49.2 8.5216899 24.6 0 0 0 0 -1235.6852

But it becomes zero, when I use the same potential for something like phonon dispersions, as shown below

('ENERGY IN PYTHON ', array([[-1409.1978995]]))

Per MPI rank memory allocation (min/avg/max) = 3.919 | 3.919 | 3.919 Mbytes

Step Temp PotEng KinEng Press Pxx Pyy

0 300 0 7.7168307 748.18204 776.72401 758.39058

I am not sure why is it happening. Can you please help me in resolving this issue?

there is not enough information here to make any helpful statement.

in general, this is a matter of debugging, so you should be able to track down how you tally the data into the global potential energy accumulators yourself. this must be an error in how your pair style is programmed.

axel.

3

Hi Axel,

Thanks for your email. I debugged this issue by using ev_tally_xyz_full() function. The energy in thermo output is matching python output. I am facing another issue related to kinetic energy calculations. The pair style is giving high jumps in temperature and pressure within first 100 steps, as shown below.

Step Temp PotEng KinEng Press Pxx Pyy c_eatoms

SCALAR IN COMPUTE PE -1409.2

SCALAR IN COMPUTE PE -1409.2

0 300 -1409.1979 7.7168307 -198314.32 -91304.419 -504347.96 -1409.1979

SCALAR IN COMPUTE PE -1741.65

100 31207.805 -1741.6546 802.75116 -1697.9005 -109611.22 24933.058 -1741.6546

I AM IN COMPUTE

I am going through the documentation of compute ke and compute press and looking for possible reasons. I will really appreciate if you can also provide insights that can help me in debugging. Below is the compute part of the pair_style I wrote. I have removed some of the code that is calling python functions for better readibility. The lammps file for the above thermo output is attached in the end. If you need any more information, please let me know and I will happy to provide.

Thanks in advance.

Regards,
Gurjot

void PairTestFinalTwo::compute(int eflag, int vflag){
setenv(“PYTHONPATH”, “.”, 0);
Py_Initialize();
PyRun_SimpleString (“import sys; sys.path.insert(0, ‘/home/gurjot/Softwares/lammps-12Dec18/src’)”);
_import_array();
std::cout<<" I AM IN COMPUTE “<<”\n";
double evdwl;
evdwl = 10.0;
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = 0;

int inum, i, ii;
int *type = atom->type;
double **x = atom->x;
double *f = atom->f;
int ilist;
double lattice;
double lp;

inum = list->inum;
ilist = list->ilist;
lattice = new double [9];
lattice[0] = double(domain->xprd);
lattice[1] = 0.0;
lattice[2] = 0.0;
lattice[3] = double(domain->xy);
lattice[4] = double(domain->yprd);
lattice[5] = 0.0;
lattice[6] = double(domain->xz);
lattice[7] = double(domain->yz);
lattice[8] = double(domain->zprd);
lp = &lattice[0];

long energyDouble;
double(*forceMatrix)[3] = new double[inum][3]; // Should work
long double(*positions)[3] = new long double[inum][3];
double *localEnergyMatrix;
localEnergyMatrix = new double [inum];
for (ii =0;ii<inum;ii++){
i = ilist[ii];
for (int jjj =0;jjj<3;jjj++){
positions[i][jjj] = x[i][jjj];
}
}

//I have deleted a bunch of code involving calls to python module
//the line below is dummy line summarizing the input and output from python code
//localEnergyMatrix is an array with each entry corresponding to energy of each atom as computed from python code
//armForce mat is armadillo matrix storing the forces for each atom in x,y,z direction
//i have debugged the python code and conversion of python to c++ and it works correctly
localEnergyMatrix,armForceMAT = getEnergyForceFromPython(positions, lattice) ;

for (int ii=0; ii<inum;ii++){
i = ilist[ii];
for (int jj =0; jj<3; jj++){
f[i][jj] = armForceMAT(i,jj);//forceMatrix[i][jj];
}
}

//POSSIBLE ERROR IN THIS PART- PRESSURE WASNT COMPUTED RIGHT
/*
eng_vdwl = 0.0;//energyDouble;
if(eflag_atom) {
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
//std::cout<<"C++ LE "<<localEnergyMatrix[ii]<<std::endl;
//std::cout<<"eatom Before "<<eatom[7]<<std::endl;
eatom[i] = localEnergyMatrix[i];
eng_vdwl = eng_vdwl + localEnergyMatrix[i];
//std::cout<<"eatom After "<<eatom[7]<<std::endl;
//if (evflag) ev_tally_xyz_full(i, localEnergyMatrix[i], 0.0, f[i][0], f[i][1], f[i][2], 1.0, 1.0, 1.0);

}
}/
if (evflag){
for (ii=0;ii<inum;ii++){
i = ilist[ii];
//had to multiply by two as energy in thermo output was coming up to be halved
ev_tally_xyz_full(i,2.0localEnergyMatrix[i],0.0, f[i][0], f[i][1], f[i][2], 1.0, 1.0, 1.0);
}
}

if (vflag_fdotr) virial_fdotr_compute();
}

Lammps file

boundary p p p

units metal

atom_style atomic

atom_modify map array

geometry

read_data data.pos

mass * 12.0107

EAM potential

#pair_style airebo 3.0

#pair_coeff * * CH.airebo C

pair_style TestFinalTwo 2.5

pair_coeff 1 1 1.0 4.2 4.2

neighbor 2. nsq

neigh_modify every 1 delay 0 check yes

---------- Define Settings ---------------------

compute eng all pe/atom

compute eatoms all pe #reduce sum c_eng

#Langevin random seed

variable dt equal 2e-3

variable r equal 57085

variable T equal 300

variable dT equal “v_dt * 100”

timestep ${dt}

initialize

velocity all create $T 28459 rot yes dist gaussian mom yes

reset_timestep 0

fixes

fix 1 all langevin $T T {dT} 73504 zero yes

fix 2 all nve

fix 3 all phonon 10 50000 3000000 map.in TestFinal0 nasr 50

output

thermo_style custom step temp pe ke press pxx pyy c_eatoms

thermo 100

#restart 1000000 restart.one restart.two

dump 1 all xyz 3000000 dump.Test.xyz

execution

run 6000000

write_restart Restart.final

4

Unless your initial condition is far from the equilibrium, this is almost always a result of an improperly implemented potential, i.e., either the potential or force is discontinuous.

5

i have no time to debug (incomplete) code for something that is looks rather opaque.

if you have sudden jumps in the kinetic energy, than that is usually due to very high forces, which means a problem in your force computation. it is very unlikely that this is due to other (commonly used) code, since those are working well with other pair styles.

6

Thanks for your valuable feedback. Working on your suggestions, I tested my pair style for few of the configurations that I used during training. Sudden jumps in KE are still there.

Parameters for this pair style were found by matching DFT values of energy and forces. During training, mean absolute error in force predictions was less than 1.0 eV/A with error in energy also of the same order. Atomic configurations used during training are strucutres from the literature and do represent the material phase.

I think I am doing something wrong when implementing it in lammps. Python code that outputs energy/forces is matching DFT values of energy and forces.

Is there a way to check discontinuity in energy/forces?

7

Yes. Implement your force and energy computation code in the “single” member function of your pair class, use pair_write to produce tabulated values of your pair potential and then plot them.