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