[gmx-users] question about calculation electrical conductivity with green-kubo relation

[Qiao Baofu]

Hi all,

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).

Is there anyone who can give some suggestions?  Thanks.


#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
#include <stdio.h>
#include <math.h>

#include "confio.h"
#include "copyrite.h"
#include "fatal.h"
#include "futil.h"
#include "gstat.h"
#include "macros.h"
#include "maths.h"
#include "physics.h"
#include "index.h"
#include "smalloc.h"
#include "statutil.h"
#include "string.h"
#include "sysstuff.h"
#include "txtdump.h"
#include "typedefs.h"
#include "vec.h"
#include "strdb.h"
#include "xvgr.h"


int main(int argc,char *argv[])
{
  static char *desc[] = {
    "This program is wroten to calculate the Electrcal Conductivity.",
    "The Green-Kubo relation is used. See the reference:",
    "J-P Hansen, I. R. McDonald: Theory for Simple Liquids, 2nd Edition,
P242."
  };
  t_pargs pa[] = {
  };
  #define NPA asize(pa)

   t_filenm  fnm[] = {
    { efTRX, "-f",  NULL,   ffREAD  },
    { efTPX, "-s",  NULL,   ffREAD  },
    { efNDX, "-n",  NULL,   ffOPTRD },
    { efXVG, "-o",  "EC",   ffWRITE }
  };
#define NFILE asize(fnm)

  FILE       *fp_EC;
  t_topology top;
  t_trxframe fr;
  bool       bTop=FALSE;
  int        gnx;
  atom_id    *index;
  char       *grpname;
  char       title[256];
  real       q,t,t0;
  int        status,i,j;
  real       *qv,*qv0;
  int        i_T_N,T_N=10,Ensemble_N;
  real       *Integrand,*Integral,EC;
  real       T_fac;

  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
  real       T=425,V=47.3,Pref;                            // Temperature,
Volume and Prefacor

  //printf("\nPlease input the Temperature(K) and Volume(nm^3),\n& How many
data points per ensemble:\n");
  //scanf("%f  %f  %g",&T,&V,&T_N);

  V=V*1.0e-27;
  Pref = 1.0/(kb*T)*(e*v_c)*(e*v_c)*t_c/V;
  printf("Prefactor = %10g\n", Pref);
  printf("\n");

  //CopyRight(stderr,argv[0]);
  parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
            NFILE,fnm,NPA,pa,asize(desc),desc,0,NULL);


bTop=read_tps_conf(ftp2fn(efTPX,NFILE,fnm),title,&top,NULL,NULL,NULL,TRUE);
  get_index(&top.atoms,ftp2fn_null(efNDX,NFILE,fnm),1,&gnx,&index,&grpname);


  sprintf(title,"Electrical Conductivity: %s",grpname);
  fp_EC=xvgropen(opt2fn("-o",NFILE,fnm),title,"t (ps)","EC (S/m)");

  snew(qv0,DIM);
  snew(qv,DIM);
  snew(Integrand,T_N);
  snew(Integral,T_N);

  read_first_frame(&status,ftp2fn(efTRX,NFILE,fnm),&fr,TRX_NEED_V);

  for(i=0; i<T_N; i++)
      Integral[i]=0;
  Ensemble_N=0;
  EC=0;

  next_ensemble:

  t0=fr.time;

  for(i=0; i<DIM; i++)
        qv0[i]=0;
  for(i=0; i<gnx; i++) {
      q = top.atoms.atom[index[i]].q;
      qv0[0] += q*fr.v[index[i]][XX];
      qv0[1] += q*fr.v[index[i]][YY];
      qv0[2] += q*fr.v[index[i]][ZZ];
  }
  i_T_N = 0;
  do {
      t=fr.time;
      for(i=0; i<DIM; i++)
           qv[i]=0;
      for(i=0; i<gnx; i++) {
          q = top.atoms.atom[index[i]].q;
          qv[0] += q*fr.v[index[i]][XX];
          qv[1] += q*fr.v[index[i]][YY];
          qv[2] += q*fr.v[index[i]][ZZ];
          }
      Integrand[i_T_N] = qv[0]*qv0[0]+qv[1]*qv0[1]+qv[2]*qv0[2];

      // T_fac is the factor between time and i_T_N
      if(i_T_N==(T_N-1))
          T_fac = (t-t0)/i_T_N;
      i_T_N++;
      if (i_T_N == T_N) {
          for (j=0; j<T_N; j++)
              Integral[j] += Integrand[j];
          Ensemble_N++;
         if (read_next_frame(status,&fr))
             goto next_ensemble;
      }
  } while (read_next_frame(status,&fr));

  close_trj(status);


  for (i=0; i<T_N; i++)
  {
      EC += Pref*Integral[i]/Ensemble_N*T_fac;
      fprintf(fp_EC, "%7.3f  %10g \n",i*T_fac,EC);
  }
  ffclose(fp_EC);
  thanx(stderr);
  return 0;
}
-- 
Sincerely yours,
**********************************************
Baofu Qiao, PhD
Frankfurt Institute for Advanced Studies
Max-von-Laue-Str. 1
60438 Frankfurt am Main, Germany TEL:+49-69-7984-7529
**********************************************
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20061030/16a4ce6a/attachment.html>