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]=delxfpair(velo[i][0]+velo[j][0]);
fv_ljc[1]=delyfpair(velo[i][1]+velo[j][1]);
fv_ljc[2]=delzfpair(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)