Hi all,<br><br>I wrote a program to calculate the electrical conductivity at zero field using the Green-Kubo relation, as follows. But I found the calculated data are strange: in the first stage, the data are very big, while in the end stage it is a negative value, also quite big. (I used the system with about 3500 atoms, 4 ns (every 100 ns are taken as one ensemble).
<br><br>Is there anyone who can give some suggestions? Thanks.<br><br><br clear="all">#ifdef HAVE_CONFIG_H<br>#include <config.h><br>#endif<br>#include <stdio.h><br>#include <math.h><br><br>#include "
confio.h"<br>#include "copyrite.h"<br>#include "fatal.h"<br>#include "futil.h"<br>#include "gstat.h"<br>#include "macros.h"<br>#include "maths.h"<br>#include "
physics.h"<br>#include "index.h"<br>#include "smalloc.h"<br>#include "statutil.h"<br>#include "string.h"<br>#include "sysstuff.h"<br>#include "txtdump.h"<br>
#include "typedefs.h"<br>#include "vec.h"<br>#include "strdb.h"<br>#include "xvgr.h"<br><br><br>int main(int argc,char *argv[])<br>{<br> static char *desc[] = {<br> "This program is wroten to calculate the Electrcal Conductivity.",
<br> "The Green-Kubo relation is used. See the reference:",<br> "J-P Hansen, I. R. McDonald: Theory for Simple Liquids, 2nd Edition, P242."<br> };<br> t_pargs pa[] = {<br> };<br> #define NPA asize(pa)
<br> <br> t_filenm fnm[] = {<br> { efTRX, "-f", NULL, ffREAD },<br> { efTPX, "-s", NULL, ffREAD }, <br> { efNDX, "-n", NULL, ffOPTRD },<br> { efXVG, "-o", "EC", ffWRITE }
<br> };<br>#define NFILE asize(fnm)<br><br> FILE *fp_EC; <br> t_topology top;<br> t_trxframe fr;<br> bool bTop=FALSE;<br> int gnx;<br> atom_id *index;<br> char *grpname;<br> char title[256];
<br> real q,t,t0;<br> int status,i,j;<br> real *qv,*qv0;<br> int i_T_N,T_N=10,Ensemble_N;<br> real *Integrand,*Integral,EC;<br> real T_fac;<br> <br> real kb=1.381e-23,e=1.6e-19
,v_c=1e3,t_c=1e-12; //v_c: the coefficient between nm/ps and m/s<br> real T=425,V=47.3,Pref; // Temperature, Volume and Prefacor <br> <br> //printf("\nPlease input the Temperature(K) and Volume(nm^3),\n& How many data points per ensemble:\n");
<br> //scanf("%f %f %g",&T,&V,&T_N);<br> <br> V=V*1.0e-27;<br> Pref = 1.0/(kb*T)*(e*v_c)*(e*v_c)*t_c/V;<br> printf("Prefactor = %10g\n", Pref); <br> printf("\n");<br> <br>
//CopyRight(stderr,argv[0]);<br> parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,<br> NFILE,fnm,NPA,pa,asize(desc),desc,0,NULL);<br> <br> bTop=read_tps_conf(ftp2fn(efTPX,NFILE,fnm),title,&top,NULL,NULL,NULL,TRUE);
<br> get_index(&top.atoms,ftp2fn_null(efNDX,NFILE,fnm),1,&gnx,&index,&grpname);<br> <br> <br> sprintf(title,"Electrical Conductivity: %s",grpname);<br> fp_EC=xvgropen(opt2fn("-o",NFILE,fnm),title,"t (ps)","EC (S/m)");
<br> <br> snew(qv0,DIM);<br> snew(qv,DIM);<br> snew(Integrand,T_N);<br> snew(Integral,T_N);<br> <br> read_first_frame(&status,ftp2fn(efTRX,NFILE,fnm),&fr,TRX_NEED_V);<br> <br> for(i=0; i<T_N; i++)<br> Integral[i]=0;
<br> Ensemble_N=0;<br> EC=0;<br> <br> next_ensemble:<br> <br> t0=fr.time;<br> <br> for(i=0; i<DIM; i++)<br> qv0[i]=0;<br> for(i=0; i<gnx; i++) {<br> q = top.atoms.atom[index[i]].q;<br> qv0[0] += q*
fr.v[index[i]][XX];<br> qv0[1] += q*fr.v[index[i]][YY];<br> qv0[2] += q*fr.v[index[i]][ZZ];<br> }<br> i_T_N = 0;<br> do {<br> t=fr.time;<br> for(i=0; i<DIM; i++)<br> qv[i]=0;<br> for(i=0; i<gnx; i++) {
<br> q = top.atoms.atom[index[i]].q;<br> qv[0] += q*fr.v[index[i]][XX];<br> qv[1] += q*fr.v[index[i]][YY];<br> qv[2] += q*fr.v[index[i]][ZZ];<br> }<br> Integrand[i_T_N] = qv[0]*qv0[0]+qv[1]*qv0[1]+qv[2]*qv0[2];
<br> <br> // T_fac is the factor between time and i_T_N<br> if(i_T_N==(T_N-1))<br> T_fac = (t-t0)/i_T_N;<br> i_T_N++;<br> if (i_T_N == T_N) {<br> for (j=0; j<T_N; j++)<br> Integral[j] += Integrand[j];
<br> Ensemble_N++;<br> if (read_next_frame(status,&fr))<br> goto next_ensemble;<br> }<br> } while (read_next_frame(status,&fr));<br> <br> close_trj(status);<br> <br><br> for (i=0; i<T_N; i++)
<br> {<br> EC += Pref*Integral[i]/Ensemble_N*T_fac;<br> fprintf(fp_EC, "%7.3f %10g \n",i*T_fac,EC);<br> }<br> ffclose(fp_EC);<br> thanx(stderr);<br> return 0;<br>}<br>-- <br>Sincerely yours,<br>**********************************************
<br>Baofu Qiao, PhD<br>Frankfurt Institute for Advanced Studies<br>Max-von-Laue-Str. 1<br>60438 Frankfurt am Main, Germany TEL:+49-69-7984-7529<br>**********************************************