Dear all:

I have gone through the compute heat flux thread by German and other users, and I found that by using the stress/atom, the original compute heat flux command is obsolete. For some reason, I need to modify the source code to compute the heat flux by myself. I started with the pair_lj_cut.cpp and calculated the heat flux (output to a separate file) due to the LJ potential and combine it with the convection term to obtain the total flux. The following is a part of the modified pair_lj_cut.cpp. The complete file is also attached. Also attached is the comparison of integrated autocorrelation function of heat flux. At the beginning, the modified lmps gives very close results to stress/atom method. However, the divergence is serious as time goes on. Is the velocities called in in pair_lj_cut.cpp is not the correct velocity to be used? or it’s just the numerical error?

//----TENGFEI MOD.

//CALCULATE VIRAL HEAT FLUX DUE TO PAIR INTERACTION

fv_ljc[0]=delx*fpair*(velo[i][0]+velo[j][0]);

fv_ljc[1]=dely*fpair*(velo[i][1]+velo[j][1]);

fv_ljc[2]=delz*fpair*(velo[i][2]+velo[j][2]);

qflux2_ljc[0]+=delx*(fv_ljc[0]+fv_ljc[1]+fv_ljc[2]);

qflux2_ljc[1]+=dely*(fv_ljc[0]+fv_ljc[1]+fv_ljc[2]);

qflux2_ljc[2]+=delz*(fv_ljc[0]+fv_ljc[1]+fv_ljc[2]);

//TENGFEI MOD.----

if (eflag) {

evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) -

offset[itype][jtype];

evdwl *= factor_lj;

}

if (evflag) ev_tally(i,j,nlocal,newton_pair,

evdwl,0.0,fpair,delx,dely,delz);

}

}

}

//----TENGFEI MOD.

//PRINT HEAT FLUX TO FILE “bond_flux”

MPI_Comm_rank(world,&me);

if (me==0){ if (nstep_ljc==1) ljc_flux=fopen(“ljc_flux”,“w”);

}

//gather all answers from other processes

MPI_Reduce(&qflux2_ljc[0], &answer[0], 3, MPI_DOUBLE, MPI_SUM, 0, world);

MPI_Comm_rank(world,&me);

//write flux to file

if (me==0) fprintf(ljc_flux,"%d %f %f %f\n",nstep_ljc,answer[0],answer[1],answer[2]);

//TENGFEI MOD.----

I hope I haven’t missed anything, more information will be provided if needed.

I appreciate any comment!

Tengfei

pair_lj_cut.cpp (22.2 KB)

in.lj.compute (879 Bytes)

in.lj.mod (1.04 KB)

in.lj.script (1.69 KB)