mdlj
 All Classes Files Functions Variables Pages
Classes | Functions
mdlj.c File Reference
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
Include dependency graph for mdlj.c:

Go to the source code of this file.

Classes

struct  dpd_t
 contains all information about the dpd thermostat More...
 
struct  vlist_t
 contains all information about the verlet list More...
 
struct  lclist_t
 contains all information about the linked cell list More...
 
struct  nblist_t
 contains all information about the neighbor list More...
 

Functions

void xyz_out (FILE *fp, double *rx, double *ry, double *rz, double *vx, double *vy, double *vz, int *ix, int *iy, int *iz, double L, int N, int z, int put_vel, int unfold)
 
int xyz_in (FILE *fp, double *rx, double *ry, double *rz, double *vx, double *vy, double *vz, int *N)
 
void save_positions (double *rx, double *ry, double *rz, int N, double *rx0, double *ry0, double *rz0)
 save current positions More...
 
double max_moved_distance2 (double *rx, double *ry, double *rz, int N, double *rx0, double *ry0, double *rz0)
 calculates the max distance a particle has moved More...
 
double per_dist2 (int i, int j, double *rx, double *ry, double *rz, double *dx, double *dy, double *dz, double L)
 calculates periodic distance between two particles More...
 
void calc_lj (int i, int j, double dx, double dy, double dz, double r2, double *fx, double *fy, double *fz, double rc2, double ecut, double *vir, double *e)
 calculates Lennard-Jones interaction between two particles More...
 
void calc_dpd (int i, int j, double dx, double dy, double dz, double r2, double *fx, double *fy, double *fz, dpd_t *dpd)
 calculates dpd force between two particles More...
 
double total_e (double *rx, double *ry, double *rz, double *fx, double *fy, double *fz, int N, double L, nblist_t *nblist, dpd_t *dpd, double rc2, double ecor, double ecut, double *vir)
 algorithm for computing forces and potential energy. algorithm for computing forces and potential energy. The virial is also computed and returned in *vir. More...
 
double drand_gaussian (double mean, double sigma)
 generate gaussian random number generate a gaussian random number using the Box-Muller method More...
 
void init (double *rx, double *ry, double *rz, double *vx, double *vy, double *vz, int *ix, int *iy, int *iz, int *N, double L, double T0, double *KE, char *icf)
 
int main (int argc, char *argv[])
 

Function Documentation

void calc_dpd ( int  i,
int  j,
double  dx,
double  dy,
double  dz,
double  r2,
double *  fx,
double *  fy,
double *  fz,
dpd_t *  dpd 
)

calculates dpd force between two particles

Parameters
[in]iindex of particle i
[in]jindex of particle j
[in]dxarray of x-coordinate of the distance vetor
[in]dyarray of y-coordinate of the distance vetor
[in]dzarray of z-coordinate of the distance vetor
[in]r2square of the distance vector
[in,out]fxarray of x-coordinate of the force of all particles
[in,out]fyarray of y-coordinate of the force of all particles
[in,out]fzarray of z-coordinate of the force of all particles
[in]dpdstructure with dpd information

Definition at line 216 of file mdlj.c.

217  {
218 
219  double v_dot_d, amplitude;
220  double tmp_x,tmp_y,tmp_z;
221  if (r2 < dpd->cutoff2) {
222  v_dot_d = ( dpd->vx[i] - dpd->vx[j] ) * dx ;
223  v_dot_d += ( dpd->vy[i] - dpd->vy[j] ) * dy ;
224  v_dot_d += ( dpd->vz[i] - dpd->vz[j] ) * dz ;
225  amplitude = dpd->prefactor / sqrt(r2) * (drand48()-0.5) - dpd->gamma * v_dot_d / r2;
226  tmp_x = amplitude * dx;
227  fx[i] += tmp_x;
228  fx[j] -= tmp_x;
229  tmp_y = amplitude * dy;
230  fy[i] += tmp_y;
231  fy[j] -= tmp_y;
232  tmp_z = amplitude * dz;
233  fz[i] += tmp_z;
234  fz[j] -= tmp_z;
235  }
236 }
dpd_t::vx
double * vx
array of x-coordinate of the velocity of all particles
Definition: mdlj.c:199
dpd_t::vz
double * vz
array of z-coordinate of the velocity of all particles
Definition: mdlj.c:201
dpd_t::prefactor
double prefactor
dpd constant = sqrt(24.0*T*gamma/dt)
Definition: mdlj.c:198
dpd_t::vy
double * vy
array of y-coordinate of the velocity of all particles
Definition: mdlj.c:200
dpd_t::gamma
double gamma
friction parameter
Definition: mdlj.c:197

Here is the caller graph for this function:

void calc_lj ( int  i,
int  j,
double  dx,
double  dy,
double  dz,
double  r2,
double *  fx,
double *  fy,
double *  fz,
double  rc2,
double  ecut,
double *  vir,
double *  e 
)

calculates Lennard-Jones interaction between two particles

Parameters
[in]iindex of particle i
[in]jindex of particle j
[in]dxarray of x-coordinate of the distance vetor
[in]dyarray of y-coordinate of the distance vetor
[in]dzarray of z-coordinate of the distance vetor
[in]r2square of the distance vector
[in,out]fxarray of x-coordinate of the force of all particles
[in,out]fyarray of y-coordinate of the force of all particles
[in,out]fzarray of z-coordinate of the force of all particles
[in]rc2square of the cutoff distance
[in]ecutenergy varlue at the cutoff
[in,out]virvariable to add the virial contribute to
[in,out]evariable to add the energy contribute to

Definition at line 175 of file mdlj.c.

177  {
178  double r6i, f;
179 
180  if (r2<rc2) {
181  r6i = 1.0/(r2*r2*r2);
182  *e += 4*(r6i*r6i - r6i) - ecut;
183  f = 48*(r6i*r6i-0.5*r6i);
184  fx[i] += dx*f/r2;
185  fx[j] -= dx*f/r2;
186  fy[i] += dy*f/r2;
187  fy[j] -= dy*f/r2;
188  fz[i] += dz*f/r2;
189  fz[j] -= dz*f/r2;
190  *vir += f;
191  }
192 }

Here is the caller graph for this function:

double drand_gaussian ( double  mean,
double  sigma 
)

generate gaussian random number generate a gaussian random number using the Box-Muller method

Parameters
[in]meanmean value of the gaussian
[in]sigmastandard deviation of the gaussian
Returns
gaussian random number

Definition at line 476 of file mdlj.c.

476  {
477  static int i=0;
478  static double x2;
479  /* http://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform */
480  double r,theta,u1,u2,x1;
481  if (i==0) {
482  i=1-i;
483  u1=1-drand48();
484  u2=1-drand48();
485  r=sqrt(-2*log(u1))*sigma;
486  theta=2*M_PI*u2;
487  x1 = r*cos(theta) + mean;
488  x2 = r*sin(theta) + mean ;
489  return(x1);
490  } else {
491  i=1-i;
492  return(x2);
493  }
494 }

Here is the caller graph for this function:

void init ( double *  rx,
double *  ry,
double *  rz,
double *  vx,
double *  vy,
double *  vz,
int *  ix,
int *  iy,
int *  iz,
int *  N,
double  L,
double  T0,
double *  KE,
char *  icf 
)

Definition at line 500 of file mdlj.c.

504  {
505  int i,iix,iiy,iiz;
506  double cmvx=0.0,cmvy=0.0,cmvz=0.0;
507  double T, fac;
508  int n3=2;
509  int vel_ok=0;
510 
511  /* If icf has a value, assume it is the name of a file containing
512  the input configuration in XYZ format */
513  if (icf) {
514  FILE * fp = fopen(icf,"r");
515  if (fp) {
516  vel_ok = xyz_in(fp,rx,ry,rz,vx,vy,vz,N);
517  }
518  else {
519  fprintf(stderr,"# error: could not read %s\n",icf);
520  exit(-1);
521  }
522  fclose(fp);
523  }
524  /* Assign particles on a cubic lattice */
525  else {
526 
527  /* Find the lowest perfect cube, n3, greater than or equal to the
528  number of particles */
529  while ((n3*n3*n3)<*N) n3++;
530 
531  iix=iiy=iiz=0;
532  /* Assign particle positions */
533  for (i=0;i<*N;i++) {
534  rx[i] = ((double)iix+0.5)*L/n3;
535  ry[i] = ((double)iiy+0.5)*L/n3;
536  rz[i] = ((double)iiz+0.5)*L/n3;
537  iix++;
538  if (iix==n3) {
539  iix=0;
540  iiy++;
541  if (iiy==n3) {
542  iiy=0;
543  iiz++;
544  }
545  }
546  }
547  }
548  /* If no velocities yet assigned, randomly pick some */
549  if (!vel_ok) {
550  for (i=0;i<*N;i++) {
551  vx[i]=drand_gaussian(0.0,1.0);
552  vy[i]=drand_gaussian(0.0,1.0);
553  vx[i]=drand_gaussian(0.0,1.0);
554  }
555  }
556  /* Take away any center-of-mass drift; compute initial KE */
557  for (i=0;i<*N;i++) {
558  cmvx+=vx[i];
559  cmvy+=vy[i];
560  cmvz+=vz[i];
561  }
562  (*KE)=0;
563  for (i=0;i<*N;i++) {
564  vx[i]-=cmvx/(*N);
565  vy[i]-=cmvy/(*N);
566  vz[i]-=cmvz/(*N);
567  (*KE)+=vx[i]*vx[i]+vy[i]*vy[i]+vz[i]*vz[i];
568  }
569  (*KE)*=0.5;
570  /* if T0 is specified, scale velocities */
571  if (T0>0.0) {
572  T=(*KE)/(*N)*2./3.;
573  fac=sqrt(T0/T);
574  (*KE)=0;
575  for (i=0;i<*N;i++) {
576  vx[i]*=fac;
577  vy[i]*=fac;
578  vz[i]*=fac;
579  (*KE)+=vx[i]*vx[i]+vy[i]*vy[i]+vz[i]*vz[i];
580  }
581  (*KE)*=0.5;
582  }
583  /* Initialize periodic boundary crossing counter arrays */
584  for (i=0;i<*N;i++) {
585  ix[i]=0;
586  iy[i]=0;
587  iz[i]=0;
588  }
589 }
drand_gaussian
double drand_gaussian(double mean, double sigma)
generate gaussian random number generate a gaussian random number using the Box-Muller method ...
Definition: mdlj.c:476
xyz_in
int xyz_in(FILE *fp, double *rx, double *ry, double *rz, double *vx, double *vy, double *vz, int *N)
Definition: mdlj.c:66

Here is the call graph for this function:

Here is the caller graph for this function:

int main ( int  argc,
char *  argv[] 
)

Definition at line 591 of file mdlj.c.

591  {
592 
593  double * rx, * ry, * rz;
594  double * vx, * vy, * vz;
595  double * fx, * fy, * fz;
596  int * ix, * iy, * iz;
597  int N=216;
598  double L=0.0;
599  double rho=pow(1./1.05,3), rc2 = 1.12246, vir, vir_old, V;
600  double PE, KE, TE, ecor, ecut, T0=1.0, TE0, T=1.0;
601  double rr3,dt=0.001, dt2;
602  int i,s;
603  int nSteps = 10, fSamp=100;
604  int use_e_corr=0;
605  int unfold = 0;
606  double langevin_gamma=0;
607  double dpd_gamma=0;
608 
609  char fn[20];
610  FILE * out;
611  char * wrt_code_str = "w";
612  char * init_cfg_file = NULL;
613  double rvl=0;
614  double rlc=0;
615 
616  unsigned long int Seed = 23410981;
617 
618  /* Here we parse the command line arguments; If
619  you add an option, document it in the usage() function! */
620  for (i=1;i<argc;i++) {
621  if (!strcmp(argv[i],"-N")) N=atoi(argv[++i]);
622  else if (!strcmp(argv[i],"-rho")) rho=atof(argv[++i]);
623  else if (!strcmp(argv[i],"-dt")) dt=atof(argv[++i]);
624  else if (!strcmp(argv[i],"-rc")) rc2=atof(argv[++i]);
625  else if (!strcmp(argv[i],"-ns")) nSteps = atoi(argv[++i]);
626  else if (!strcmp(argv[i],"-T0")) T0=atof(argv[++i]);
627  else if (!strcmp(argv[i],"-T")) T=atof(argv[++i]);
628  else if (!strcmp(argv[i],"-sd")) langevin_gamma=atof(argv[++i]);
629  else if (!strcmp(argv[i],"-dpd")) dpd_gamma=atof(argv[++i]);
630  else if (!strcmp(argv[i],"-fs")) fSamp=atoi(argv[++i]);
631  else if (!strcmp(argv[i],"-sf")) wrt_code_str = argv[++i];
632  else if (!strcmp(argv[i],"-icf")) init_cfg_file = argv[++i];
633  else if (!strcmp(argv[i],"-ecorr")) use_e_corr = 1;
634  else if (!strcmp(argv[i],"-seed")) Seed = (unsigned long)atoi(argv[++i]);
635  else if (!strcmp(argv[i],"-rvl")) rvl = atof(argv[++i]);
636  else if (!strcmp(argv[i],"-rlc")) rlc = atof(argv[++i]);
637  else if (!strcmp(argv[i],"-uf")) unfold = 1;
638  else if ((!strcmp(argv[i],"-h"))||(!strcmp(argv[i],"--help"))) {
639  fprintf(stdout,"mdlj usage:\n");
640  fprintf(stdout,"mdlj [options]\n\n");
641  fprintf(stdout,"Options:\n");
642  fprintf(stdout,"\t -N [integer]\t\tNumber of particles (default %i)\n",N);
643  fprintf(stdout,"\t -rho [real]\t\tNumber density (default %f)\n",rho);
644  fprintf(stdout,"\t -dt [real]\t\tTime step (default %f)\n",dt);
645  fprintf(stdout,"\t -rc [real]\t\tCutoff radius (default %f)\n",rc2);
646  fprintf(stdout,"\t -ns [real]\t\tNumber of integration steps (default %i)\n",nSteps);
647  fprintf(stdout,"\t -T0 [real]\t\tInitial temperature (default %f)\n",T0);
648  fprintf(stdout,"\t -sd [real]\t\tUse Langevin thermostat with frition constant 0=off (default %f)\n",langevin_gamma);
649  fprintf(stdout,"\t -dpd [real]\t\tUse DPD thermostat with frition constant 0=off (default %f)\n",dpd_gamma);
650  fprintf(stdout,"\t -T [real]\t\tTemperature of the Langevin thermostat(default %f)\n",T);
651  fprintf(stdout,"\t -fs [integer]\t\tSample frequency (default %i)\n",fSamp);
652  fprintf(stdout,"\t -sf [a|w]\t\tAppend or write config output file (default %s)\n",wrt_code_str);
653  fprintf(stdout,"\t -icf [string]\t\tInitial configuration file (default %s)\n",init_cfg_file);
654  fprintf(stdout,"\t -seed [integer]\tRandom number generator seed (default %li)\n",Seed);
655  fprintf(stdout,"\t -rvl [real]\t\tUse neighbor lists with this cutoff (default %f)\n",rvl);
656  fprintf(stdout,"\t -rlc [read] \t\tUse linked cells with this size (default %f)\n",rlc);
657  fprintf(stdout,"\t -uf \t\tPrint unfolded coordinates in output files (default %i)\n",unfold);
658  fprintf(stdout,"\t -h \t\tPrint this info\n");
659  exit(0);
660  }
661  else {
662  fprintf(stderr,"Error: Command-line argument '%s' not recognized.\n",
663  argv[i]);
664  exit(-1);
665  }
666  }
667 
668  /* Compute the side-length */
669  L = pow((V=N/rho),0.3333333);
670 
671  /* Compute the tail-corrections; assumes sigma and epsilon are both 1 */
672  rr3 = 1.0/(rc2*rc2*rc2);
673  ecor = use_e_corr?8*M_PI*rho*(rr3*rr3*rr3/9.0-rr3/3.0):0.0;
674  ecut = 4*(rr3*rr3*rr3*rr3-rr3*rr3);
675 
676  /* Compute the *squared* cutoff, reusing the variable rc2 */
677  rc2*=rc2;
678 
679  /* compute the squared time step */
680  dt2=dt*dt;
681 
682  /* Output some initial information */
683  fprintf(stdout,"# NVE MD Simulation of a Lennard-Jones fluid\n");
684  fprintf(stdout,"# L = %.5lf; rho = %.5lf; N = %i; rc = %.5lf\n",
685  L,rho,N,sqrt(rc2));
686  fprintf(stdout,"# nSteps %i, seed %li, dt %.5lf\n", nSteps,Seed,dt);
687 
688  nblist_t *nblist=NULL;
689 
690  if (( rvl > 0 ) || ( rlc > 0 )){
691  nblist=(nblist_t *)malloc(sizeof(nblist_t));
692  }
693 
694  /* verlet list stuff */
695  if (rvl > 0) {
696  nblist->vl=(vlist_t *)malloc(sizeof(vlist_t));
697  nblist->vl->cutoff = rvl;
698  nblist->vl->cutoff2 = rvl*rvl;
699  if (nblist->vl->cutoff2 < rc2) {
700  fprintf(stderr,"-rvl must be bigger than -rc\n");
701  exit(1);
702  }
703  fprintf(stdout,"# using Verlet lists with cutoff %.5f\n",nblist->vl->cutoff);
704  /* assume 25 neighbors and realloc later if needed */
705  nblist->vl->skin2 = pow(rvl-sqrt(rc2),2);
706  nblist->vl->size= 25*N;
707  nblist->vl->neighbors = (int *)malloc(nblist->vl->size*sizeof(int));
708  nblist->vl->head = (int *)malloc(N*sizeof(int));
709  /* always update nblist one time */
710  nblist->vl->update=1;
711  /* old position saver*/
712  nblist->vl->rx0 = (double*)malloc(N*sizeof(double));
713  nblist->vl->ry0 = (double*)malloc(N*sizeof(double));
714  nblist->vl->rz0 = (double*)malloc(N*sizeof(double));
715  }
716 
717  if ( rlc > 0 ) {
718  nblist->lc=(lclist_t *)malloc(sizeof(lclist_t));
719  nblist->lc->cutoff = rlc;
720  nblist->lc->ncells=(int)floor(L/rlc);
721  nblist->lc->ncells2=nblist->lc->ncells * nblist->lc->ncells;
722  nblist->lc->ncells3=nblist->lc->ncells * nblist->lc->ncells2;
723  nblist->lc->left=nblist->lc->right=1;
724  if (nblist->lc->ncells==1){
725  nblist->lc->left=nblist->lc->right=0;
726  }
727  else if (nblist->lc->ncells==2){
728  nblist->lc->left=0;
729  }
730  nblist->lc->lcell = L/nblist->lc->ncells;
731  nblist->lc->skin2 = pow(nblist->lc->lcell-sqrt(rc2),2);
732  if (rlc*rlc < rc2) {
733  fprintf(stderr,"-rlc must be bigger than -rc\n");
734  exit(1);
735  }
736  fprintf(stdout,"# using Linked cell lists with grid %i x %i x %i (cufoff %f)\n",nblist->lc->ncells,nblist->lc->ncells,nblist->lc->ncells,nblist->lc->lcell);
737  nblist->lc->neighbors = (int *)malloc(N*sizeof(int));
738  nblist->lc->head = (int *)malloc(nblist->lc->ncells3*sizeof(int));
739  /* always update nblist one time */
740  nblist->lc->update=1;
741  /* old position saver*/
742  nblist->lc->rx0 = (double*)malloc(N*sizeof(double));
743  nblist->lc->ry0 = (double*)malloc(N*sizeof(double));
744  nblist->lc->rz0 = (double*)malloc(N*sizeof(double));
745  if(nblist->vl) {
746  nblist->vl->id_head = (int *)malloc(N*sizeof(int));
747  /* we need to store an end head for the Nth particle, too */
748  nblist->vl->head = (int *)realloc(nblist->vl->head,(N+1)*sizeof(int));
749  }
750  }
751 
752  if(!nblist){
753  fprintf(stdout,"# no neighbor lists\n");
754  }
755 
756  double langevin_const = 0;
757  if (langevin_gamma > 0) {
758  fprintf(stdout,"# using Langevin thermostat for T=%f with gamma=%f\n",T,langevin_gamma);
759  langevin_const = sqrt(24.0*T*langevin_gamma/dt);
760  if (dpd_gamma > 0 ) {
761  fprintf(stderr,"# error: dpd and langevin thermostat cannot be enable at the same time\n");
762  exit(-1);
763  }
764  }
765 
766  /* Seed the random number generator */
767  srand48(Seed);
768 
769  /* Allocate the position arrays */
770  rx = (double*)malloc(N*sizeof(double));
771  ry = (double*)malloc(N*sizeof(double));
772  rz = (double*)malloc(N*sizeof(double));
773 
774  /* Allocate the boundary crossing counter arrays */
775  ix = (int*)malloc(N*sizeof(int));
776  iy = (int*)malloc(N*sizeof(int));
777  iz = (int*)malloc(N*sizeof(int));
778 
779  /* Allocate the velocity arrays */
780  vx = (double*)malloc(N*sizeof(double));
781  vy = (double*)malloc(N*sizeof(double));
782  vz = (double*)malloc(N*sizeof(double));
783 
784  /* Allocate the force arrays */
785  fx = (double*)malloc(N*sizeof(double));
786  fy = (double*)malloc(N*sizeof(double));
787  fz = (double*)malloc(N*sizeof(double));
788 
789  /* dpd thermostat */
790  dpd_t *dpd = NULL;
791  if(dpd_gamma > 0) {
792  dpd = (dpd_t *)malloc(sizeof(dpd_t));
793  dpd->gamma = dpd_gamma;
794  dpd->cutoff2 = rc2; /*could be different */
795  dpd->prefactor = sqrt(24.0*T*dpd_gamma/dt);
796  dpd->vx = vx;
797  dpd->vy = vy;
798  dpd->vz = vz;
799  }
800 
801  /* Generate initial positions on a cubic grid,
802  and measure initial energy */
803  init(rx,ry,rz,vx,vy,vz,ix,iy,iz,&N,L,T0,&KE,init_cfg_file);
804  sprintf(fn,"%i.xyz",0);
805  out=fopen(fn,"w");
806  xyz_out(out,rx,ry,rz,vx,vy,vz,ix,iy,iz,L,N,16,1,unfold);
807  fclose(out);
808 
809  if (nblist){
810  if (nblist->vl) {
811  save_positions(rx,ry,rz,N,nblist->vl->rx0,nblist->vl->ry0,nblist->vl->rz0);
812  }
813  if (nblist->lc) {
814  save_positions(rx,ry,rz,N,nblist->lc->rx0,nblist->lc->ry0,nblist->lc->rz0);
815  }
816  }
817 
818  /* DPD thermostat */
819  /* in the initial md step the thermostat only acts over "half" dt */
820  if (dpd) {
821  dpd->gamma = dpd_gamma * sqrt(3);
822  PE = total_e(rx,ry,rz,fx,fy,fz,N,L,nblist,dpd,rc2,ecor,ecut,&vir_old);
823  dpd->gamma = dpd_gamma;
824  } else {
825  PE = total_e(rx,ry,rz,fx,fy,fz,N,L,nblist,dpd,rc2,ecor,ecut,&vir_old);
826  }
827  TE0=PE+KE;
828 
829  /* Langevin thermostat */
830  /* in the initial md step the thermostat only acts over "half" dt */
831  if (langevin_gamma > 0) {
832  for (i=0;i<N;i++) {
833  fx[i]+= -langevin_gamma*vx[i] + langevin_const*sqrt(3)*(drand48()-0.5);
834  fy[i]+= -langevin_gamma*vy[i] + langevin_const*sqrt(3)*(drand48()-0.5);
835  fz[i]+= -langevin_gamma*vz[i] + langevin_const*sqrt(3)*(drand48()-0.5);
836  }
837  }
838 
839  fprintf(stdout,"# step time PotE KinE TotE Tdrift T P\n");
840 
841  for (s=0;s<nSteps;s++) {
842 
843  /* First integration half-step */
844  for (i=0;i<N;i++) {
845  rx[i]+=vx[i]*dt+0.5*dt2*fx[i];
846  ry[i]+=vy[i]*dt+0.5*dt2*fy[i];
847  rz[i]+=vz[i]*dt+0.5*dt2*fz[i];
848  vx[i]+=0.5*dt*fx[i];
849  vy[i]+=0.5*dt*fy[i];
850  vz[i]+=0.5*dt*fz[i];
851  /* Apply periodic boundary conditions */
852  if (rx[i]<0.0) { rx[i]+=L; ix[i]--; }
853  if (rx[i]>L) { rx[i]-=L; ix[i]++; }
854  if (ry[i]<0.0) { ry[i]+=L; iy[i]--; }
855  if (ry[i]>L) { ry[i]-=L; iy[i]++; }
856  if (rz[i]<0.0) { rz[i]+=L; iz[i]--; }
857  if (rz[i]>L) { rz[i]-=L; iz[i]++; }
858  }
859 
860  if (nblist) {
861  /* auto nblist update only */
862  /*2* dx > skin */
863  if ((nblist->lc)&&(nblist->lc->update==0)&&
864  (4*max_moved_distance2(rx,ry,rz,N,nblist->lc->rx0,nblist->lc->ry0,nblist->lc->rz0) > nblist->lc->skin2)) {
865  nblist->lc->update=1;
866  save_positions(rx,ry,rz,N,nblist->lc->rx0,nblist->lc->ry0,nblist->lc->rz0);
867  /* trigger verlet list update too */
868  if (nblist->vl){
869  nblist->vl->update=1;
870  save_positions(rx,ry,rz,N,nblist->vl->rx0,nblist->vl->ry0,nblist->vl->rz0);
871  }
872  }
873  if ((nblist->vl)&&(nblist->vl->update==0)&&
874  (4*max_moved_distance2(rx,ry,rz,N,nblist->vl->rx0,nblist->vl->ry0,nblist->vl->rz0) > nblist->vl->skin2)) {
875  nblist->vl->update=1;
876  save_positions(rx,ry,rz,N,nblist->vl->rx0,nblist->vl->ry0,nblist->vl->rz0);
877  }
878  }
879 
880  /* Calculate forces */
881  PE = total_e(rx,ry,rz,fx,fy,fz,N,L,nblist,dpd,rc2,ecor,ecut,&vir);
882 
883  /* Langevin thermostat */
884  if (langevin_gamma > 0) {
885  for (i=0;i<N;i++) {
886  fx[i]+= -langevin_gamma*vx[i] + langevin_const*(drand48()-0.5);
887  fy[i]+= -langevin_gamma*vy[i] + langevin_const*(drand48()-0.5);
888  fz[i]+= -langevin_gamma*vz[i] + langevin_const*(drand48()-0.5);
889  }
890  }
891 
892  /* Second integration half-step */
893  KE = 0.0;
894  for (i=0;i<N;i++) {
895  vx[i]+=0.5*dt*fx[i];
896  vy[i]+=0.5*dt*fy[i];
897  vz[i]+=0.5*dt*fz[i];
898  KE+=vx[i]*vx[i]+vy[i]*vy[i]+vz[i]*vz[i];
899  }
900  KE*=0.5;
901  TE=PE+KE;
902  fprintf(stdout,"%i %.5lf %.5lf %.5lf %.5lf %.5le %.5lf %.5lf\n",
903  s,s*dt,PE,KE,TE,(TE-TE0)/TE0,KE*2/3./N,rho*KE*2./3./N+vir/3.0/V);
904  if (!(s%fSamp)) {
905  sprintf(fn,"%i.xyz",!strcmp(wrt_code_str,"a")?0:s);
906  out=fopen(fn,wrt_code_str);
907  xyz_out(out,rx,ry,rz,vx,vy,vz,ix,iy,iz,L,N,16,1,unfold);
908  fclose(out);
909  }
910  }
911  if(nblist) {
912  if(nblist->lc){
913  free(nblist->lc->neighbors);
914  free(nblist->lc->head);
915  free(nblist->lc->rx0);
916  free(nblist->lc->ry0);
917  free(nblist->lc->rz0);
918  free(nblist->lc);
919  if(nblist->vl){
920  free(nblist->vl->id_head);
921  }
922  }
923  if(nblist->vl){
924  free(nblist->vl->neighbors);
925  free(nblist->vl->head);
926  free(nblist->vl->rx0);
927  free(nblist->vl->ry0);
928  free(nblist->vl->rz0);
929  free(nblist->vl);
930  }
931  free(nblist);
932  }
933  free(rx);
934  free(ry);
935  free(rz);
936  free(vx);
937  free(vy);
938  free(vz);
939  free(ix);
940  free(iy);
941  free(iz);
942  free(fx);
943  free(fy);
944  free(fz);
945  return(0);
946 }
lclist_t::ncells
int ncells
number of cell in each direction
Definition: mdlj.c:257
lclist_t::ncells2
int ncells2
square of number of cell in each direction
Definition: mdlj.c:258
vlist_t::cutoff
double cutoff
verlet list cutoff
Definition: mdlj.c:244
vlist_t::rx0
double * rx0
array of x-coordinate of the old position of all particles
Definition: mdlj.c:247
vlist_t::update
int update
if verlet list should be updated
Definition: mdlj.c:243
lclist_t::update
int update
if verlet list should be updated
Definition: mdlj.c:256
lclist_t::right
int right
if there is a neighbor on the right
Definition: mdlj.c:261
dpd_t::vx
double * vx
array of x-coordinate of the velocity of all particles
Definition: mdlj.c:199
init
void init(double *rx, double *ry, double *rz, double *vx, double *vy, double *vz, int *ix, int *iy, int *iz, int *N, double L, double T0, double *KE, char *icf)
Definition: mdlj.c:500
lclist_t::ncells3
int ncells3
total of number of cells (ncells**3)
Definition: mdlj.c:259
lclist_t
contains all information about the linked cell list
Definition: mdlj.c:253
xyz_out
void xyz_out(FILE *fp, double *rx, double *ry, double *rz, double *vx, double *vy, double *vz, int *ix, int *iy, int *iz, double L, int N, int z, int put_vel, int unfold)
Definition: mdlj.c:47
vlist_t
contains all information about the verlet list
Definition: mdlj.c:238
vlist_t::head
int * head
point to first neighboring particles of the ith particle
Definition: mdlj.c:239
lclist_t::neighbors
int * neighbors
list of all neighboring paricles
Definition: mdlj.c:255
max_moved_distance2
double max_moved_distance2(double *rx, double *ry, double *rz, int N, double *rx0, double *ry0, double *rz0)
calculates the max distance a particle has moved
Definition: mdlj.c:110
vlist_t::neighbors
int * neighbors
list of all neighboring paricles
Definition: mdlj.c:241
lclist_t::lcell
double lcell
length of a cell (effective cutoff)
Definition: mdlj.c:262
lclist_t::left
int left
if there is a neighbor on the left
Definition: mdlj.c:260
lclist_t::ry0
double * ry0
array of y-coordinate of the old position of all particles
Definition: mdlj.c:266
vlist_t::skin2
double skin2
square of the verlet list skin (list cutoff - lj cutoff)
Definition: mdlj.c:246
lclist_t::head
int * head
point to first neighboring particles of the ith cell
Definition: mdlj.c:254
lclist_t::rx0
double * rx0
array of x-coordinate of the old position of all particles
Definition: mdlj.c:265
vlist_t::id_head
int * id_head
particle id of particle in head needed in combination with linked cell
Definition: mdlj.c:240
nblist_t::lc
lclist_t * lc
pointer to linked cell list struct
Definition: mdlj.c:273
vlist_t::ry0
double * ry0
array of y-coordinate of the old position of all particles
Definition: mdlj.c:248
lclist_t::skin2
double skin2
square of the linked cell list skin (list cutoff - lj cutoff)
Definition: mdlj.c:264
nblist_t::vl
vlist_t * vl
pointer to verlet list struct
Definition: mdlj.c:272
vlist_t::cutoff2
double cutoff2
square of the verlet list cutoff
Definition: mdlj.c:245
vlist_t::size
int size
size of neighbors array
Definition: mdlj.c:242
nblist_t
contains all information about the neighbor list
Definition: mdlj.c:271
dpd_t::vz
double * vz
array of z-coordinate of the velocity of all particles
Definition: mdlj.c:201
dpd_t
contains all information about the dpd thermostat
Definition: mdlj.c:195
dpd_t::prefactor
double prefactor
dpd constant = sqrt(24.0*T*gamma/dt)
Definition: mdlj.c:198
total_e
double total_e(double *rx, double *ry, double *rz, double *fx, double *fy, double *fz, int N, double L, nblist_t *nblist, dpd_t *dpd, double rc2, double ecor, double ecut, double *vir)
algorithm for computing forces and potential energy. algorithm for computing forces and potential ene...
Definition: mdlj.c:280
lclist_t::rz0
double * rz0
array of z-coordinate of the old position of all particles
Definition: mdlj.c:267
save_positions
void save_positions(double *rx, double *ry, double *rz, int N, double *rx0, double *ry0, double *rz0)
save current positions
Definition: mdlj.c:90
dpd_t::vy
double * vy
array of y-coordinate of the velocity of all particles
Definition: mdlj.c:200
vlist_t::rz0
double * rz0
array of z-coordinate of the old position of all particles
Definition: mdlj.c:249
lclist_t::cutoff
double cutoff
list cutoff, round up from user input
Definition: mdlj.c:263
dpd_t::cutoff2
double cutoff2
square of the cutoff
Definition: mdlj.c:196
dpd_t::gamma
double gamma
friction parameter
Definition: mdlj.c:197

Here is the call graph for this function:

double max_moved_distance2 ( double *  rx,
double *  ry,
double *  rz,
int  N,
double *  rx0,
double *  ry0,
double *  rz0 
)

calculates the max distance a particle has moved

Parameters
[in]rxarray of x-coordinate of the position of all particles
[in]ryarray of y-coordinate of the position of all particles
[in]rzarray of z-coordinate of the position of all particles
[in]Nnumber of particles
[in]rx0array of x-coordinate of the position of all particles
[in]ry0array of y-coordinate of the position of all particles
[in]rz0array of z-coordinate of the position of all particles
Returns
square of the max distance moved

Definition at line 110 of file mdlj.c.

111  {
112  double max2=0;
113  int i;
114  for(i=0;i<N;i++) {
115  double dx2= (rx[i]-rx0[i])*(rx[i]-rx0[i]) \
116  +(ry[i]-ry0[i])*(ry[i]-ry0[i]) \
117  +(rz[i]-rz0[i])*(rz[i]-rz0[i]);
118  if (dx2>max2) {
119  max2=dx2;
120  }
121  }
122  return max2;
123 }

Here is the caller graph for this function:

double per_dist2 ( int  i,
int  j,
double *  rx,
double *  ry,
double *  rz,
double *  dx,
double *  dy,
double *  dz,
double  L 
)

calculates periodic distance between two particles

Parameters
[in]iindex of particle i
[in]jindex of particle j
[in]rxarray of x-coordinate of the position of all particles
[in]ryarray of y-coordinate of the position of all particles
[in]rzarray of z-coordinate of the position of all particles
[out]dxdifference in x-direction
[out]dydifference in y-direction
[out]dzdifference in z-direction
[in]Lbox length
Returns
square shorted distance of all periodic images

Definition at line 137 of file mdlj.c.

139  {
140  double r2,hL=L/2.0;
141 
142  /* Periodic boundary conditions: Apply the minimum image
143  convention; note that this is *not* used to truncate the
144  potential as long as there an explicit cutoff. */
145 
146  *dx = (rx[i]-rx[j]);
147  *dy = (ry[i]-ry[j]);
148  *dz = (rz[i]-rz[j]);
149 
150  if (*dx>hL) *dx-=L;
151  else if (*dx<-hL) *dx+=L;
152  if (*dy>hL) *dy-=L;
153  else if (*dy<-hL) *dy+=L;
154  if (*dz>hL) *dz-=L;
155  else if (*dz<-hL) *dz+=L;
156  r2 = *dx* *dx + *dy* *dy + *dz* *dz;
157  return r2;
158 }

Here is the caller graph for this function:

void save_positions ( double *  rx,
double *  ry,
double *  rz,
int  N,
double *  rx0,
double *  ry0,
double *  rz0 
)

save current positions

Parameters
[in]rxarray of x-coordinate of the position of all particles
[in]ryarray of y-coordinate of the position of all particles
[in]rzarray of z-coordinate of the position of all particles
[in]Nnumber of particles
[out]rx0array of x-coordinate of the position of all particles
[out]ry0array of y-coordinate of the position of all particles
[out]rz0array of z-coordinate of the position of all particles

Definition at line 90 of file mdlj.c.

91  {
92  int i;
93  for (i=0;i<N;i++) {
94  rx0[i]=rx[i];
95  ry0[i]=ry[i];
96  rz0[i]=rz[i];
97  }
98 }

Here is the caller graph for this function:

double total_e ( double *  rx,
double *  ry,
double *  rz,
double *  fx,
double *  fy,
double *  fz,
int  N,
double  L,
nblist_t *  nblist,
dpd_t *  dpd,
double  rc2,
double  ecor,
double  ecut,
double *  vir 
)

algorithm for computing forces and potential energy. algorithm for computing forces and potential energy. The virial is also computed and returned in *vir.

Definition at line 280 of file mdlj.c.

283  {
284  int i,j,k,l;
285  int cx,cy,cz,nbx,nby,nbz,c;
286  double e = 0.0;
287 
288  double dx,dy,dz,r2;
289 
290  *vir=0.0;
291 
292  /* Zero the forces */
293  for (i=0;i<N;i++) {
294  fx[i]=fy[i]=fz[i]=0.0;
295  }
296 
297  if (!nblist) {
298  /* do N square calculation */
299  for (i=0;i<(N-1);i++) {
300  for (j=i+1;j<N;j++) {
301  r2=per_dist2(i,j,rx,ry,rz,&dx,&dy,&dz,L);
302  calc_lj(i,j,dx,dy,dz,r2,fx,fy,fz,rc2,ecut,vir,&e);
303  if (dpd){
304  calc_dpd(i,j,dx,dy,dz,r2,fx,fy,fz,dpd);
305  }
306  }
307  }
308  return e+N*ecor;
309  }
310 
311  /* rebuild cell list */
312  if (nblist->lc && nblist->lc->update) {
313  for(c=0;c<nblist->lc->ncells3;c++){
314  nblist->lc->head[c]=-1;
315  }
316  /* this is an array implementation of a linked list*/
317  for(i=0;i<N;i++) {
318  cx=(int)floor(rx[i]/nblist->lc->lcell);
319  cy=(int)floor(ry[i]/nblist->lc->lcell);
320  cz=(int)floor(rz[i]/nblist->lc->lcell);
321  c=cx+cy*nblist->lc->ncells+cz*nblist->lc->ncells2;
322  nblist->lc->neighbors[i] = nblist->lc->head[c];
323  nblist->lc->head[c] = i;
324  }
325  nblist->lc->update=0;
326  }
327 
328  /* no verlet list, simple update */
329  if(!nblist->vl){
330  /* loop of cells in 3dim */
331  for(cx=0;cx<nblist->lc->ncells;cx++){
332  for(cy=0;cy<nblist->lc->ncells;cy++){
333  for(cz=0;cz<nblist->lc->ncells;cz++){
334  c=cx+cy*nblist->lc->ncells+cz*nblist->lc->ncells2;
335  /* loop over all neighbor cell of c */
336  for (nbx=cx-nblist->lc->left;nbx<=cx+nblist->lc->right;nbx++){
337  for (nby=cy-nblist->lc->left;nby<=cy+nblist->lc->right;nby++){
338  for (nbz=cz-nblist->lc->left;nbz<=cz+nblist->lc->right;nbz++){
339  /* calc cell number with respect to periodicity */
340  int nb=((nbx+nblist->lc->ncells)%nblist->lc->ncells)
341  +((nby+nblist->lc->ncells)%nblist->lc->ncells)*nblist->lc->ncells
342  +((nbz+nblist->lc->ncells)%nblist->lc->ncells)*nblist->lc->ncells2;
343  for (i=nblist->lc->head[c];i!=-1;i = nblist->lc->neighbors[i]) {
344  for (j=nblist->lc->head[nb];j != -1;j = nblist->lc->neighbors[j]) {
345  if ( i < j ) {
346  r2=per_dist2(i,j,rx,ry,rz,&dx,&dy,&dz,L);
347  calc_lj(i,j,dx,dy,dz,r2,fx,fy,fz,rc2,ecut,vir,&e);
348  if (dpd){
349  calc_dpd(i,j,dx,dy,dz,r2,fx,fy,fz,dpd);
350  }
351  }
352  }
353  }
354  }
355  }
356  }
357  }
358  }
359  }
360  return e+N*ecor;
361  }
362 
363  /* verlet list update from cell lists*/
364  if (nblist->lc && nblist->vl->update) {
365  k=0;
366  l=0;
367  /* loop of cells in 3dim */
368  for(cx=0;cx<nblist->lc->ncells;cx++){
369  for(cy=0;cy<nblist->lc->ncells;cy++){
370  for(cz=0;cz<nblist->lc->ncells;cz++){
371  c=cx+cy*nblist->lc->ncells+cz*nblist->lc->ncells2;
372  /* note different loop order than above*/
373  for (i=nblist->lc->head[c];i!=-1;i=nblist->lc->neighbors[i]) {
374  nblist->vl->head[l]=k;
375  nblist->vl->id_head[l++]=i;
376  /* loop over all neighbor cell of c */
377  for (nbx=cx-nblist->lc->left;nbx<=cx+nblist->lc->right;nbx++){
378  for (nby=cy-nblist->lc->left;nby<=cy+nblist->lc->right;nby++){
379  for (nbz=cz-nblist->lc->left;nbz<=cz+nblist->lc->right;nbz++){
380  /* calc cell number with respect to periodicity */
381  int nb=((nbx+nblist->lc->ncells)%nblist->lc->ncells)
382  +((nby+nblist->lc->ncells)%nblist->lc->ncells)*nblist->lc->ncells
383  +((nbz+nblist->lc->ncells)%nblist->lc->ncells)*nblist->lc->ncells2;
384  for (j=nblist->lc->head[nb];j != -1;j = nblist->lc->neighbors[j]) {
385  if ( i < j ) {
386  r2=per_dist2(i,j,rx,ry,rz,&dx,&dy,&dz,L);
387  if (r2 < nblist->vl->cutoff2) {
388  if(k == nblist->vl->size ) {
389  /* double the size */
390  nblist->vl->size*=2;
391  nblist->vl->neighbors=(int *)realloc(nblist->vl->neighbors,nblist->vl->size*sizeof(int));
392  }
393  nblist->vl->neighbors[k++]=j;
394  calc_lj(i,j,dx,dy,dz,r2,fx,fy,fz,rc2,ecut,vir,&e);
395  if (dpd){
396  calc_dpd(i,j,dx,dy,dz,r2,fx,fy,fz,dpd);
397  }
398  }
399  }
400  }
401  }
402  }
403  }
404  }
405  }
406  }
407  }
408  nblist->vl->head[l]=k;
409  nblist->vl->update=0;
410  return e+N*ecor;
411  }
412 
413  /* verlet list update without cell list */
414  if (nblist->vl->update) {
415  /* update neighbor list and calc lj */
416  k=0;
417  for (i=0;i<(N-1);i++) {
418  nblist->vl->head[i]=k;
419  for (j=i+1;j<N;j++) {
420  r2=per_dist2(i,j,rx,ry,rz,&dx,&dy,&dz,L);
421  if (r2 < nblist->vl->cutoff2) {
422  if(k == nblist->vl->size ) {
423  /* double the size */
424  nblist->vl->size*=2;
425  nblist->vl->neighbors=(int *)realloc(nblist->vl->neighbors,nblist->vl->size*sizeof(int));
426  }
427  nblist->vl->neighbors[k++]=j;
428  calc_lj(i,j,dx,dy,dz,r2,fx,fy,fz,rc2,ecut,vir,&e);
429  if (dpd){
430  calc_dpd(i,j,dx,dy,dz,r2,fx,fy,fz,dpd);
431  }
432  }
433  }
434  }
435  nblist->vl->head[N-1]=k;
436  nblist->vl->update=0;
437  return e+N*ecor;
438  }
439 
440  /* use existing neighbor list with linked id mapping */
441  if (nblist->vl->id_head){
442  for (l=0;l<N;l++) {
443  i=nblist->vl->id_head[l];
444  for (k=nblist->vl->head[l];k<nblist->vl->head[l+1];k++) {
445  j=nblist->vl->neighbors[k];
446  r2=per_dist2(i,j,rx,ry,rz,&dx,&dy,&dz,L);
447  calc_lj(i,j,dx,dy,dz,r2,fx,fy,fz,rc2,ecut,vir,&e);
448  if (dpd){
449  calc_dpd(i,j,dx,dy,dz,r2,fx,fy,fz,dpd);
450  }
451  }
452  }
453  } else {
454  /* use existing neighbor list verlit list only */
455  for (i=0;i<(N-1);i++) {
456  for (k=nblist->vl->head[i];k<nblist->vl->head[i+1];k++) {
457  j=nblist->vl->neighbors[k];
458  r2=per_dist2(i,j,rx,ry,rz,&dx,&dy,&dz,L);
459  calc_lj(i,j,dx,dy,dz,r2,fx,fy,fz,rc2,ecut,vir,&e);
460  if (dpd){
461  calc_dpd(i,j,dx,dy,dz,r2,fx,fy,fz,dpd);
462  }
463  }
464  }
465  }
466  return e+N*ecor;
467 }
calc_dpd
void calc_dpd(int i, int j, double dx, double dy, double dz, double r2, double *fx, double *fy, double *fz, dpd_t *dpd)
calculates dpd force between two particles
Definition: mdlj.c:216
lclist_t::ncells
int ncells
number of cell in each direction
Definition: mdlj.c:257
lclist_t::ncells2
int ncells2
square of number of cell in each direction
Definition: mdlj.c:258
vlist_t::update
int update
if verlet list should be updated
Definition: mdlj.c:243
lclist_t::update
int update
if verlet list should be updated
Definition: mdlj.c:256
lclist_t::ncells3
int ncells3
total of number of cells (ncells**3)
Definition: mdlj.c:259
per_dist2
double per_dist2(int i, int j, double *rx, double *ry, double *rz, double *dx, double *dy, double *dz, double L)
calculates periodic distance between two particles
Definition: mdlj.c:137
vlist_t::head
int * head
point to first neighboring particles of the ith particle
Definition: mdlj.c:239
lclist_t::neighbors
int * neighbors
list of all neighboring paricles
Definition: mdlj.c:255
vlist_t::neighbors
int * neighbors
list of all neighboring paricles
Definition: mdlj.c:241
lclist_t::lcell
double lcell
length of a cell (effective cutoff)
Definition: mdlj.c:262
lclist_t::left
int left
if there is a neighbor on the left
Definition: mdlj.c:260
lclist_t::head
int * head
point to first neighboring particles of the ith cell
Definition: mdlj.c:254
vlist_t::id_head
int * id_head
particle id of particle in head needed in combination with linked cell
Definition: mdlj.c:240
nblist_t::lc
lclist_t * lc
pointer to linked cell list struct
Definition: mdlj.c:273
nblist_t::vl
vlist_t * vl
pointer to verlet list struct
Definition: mdlj.c:272
vlist_t::size
int size
size of neighbors array
Definition: mdlj.c:242
calc_lj
void calc_lj(int i, int j, double dx, double dy, double dz, double r2, double *fx, double *fy, double *fz, double rc2, double ecut, double *vir, double *e)
calculates Lennard-Jones interaction between two particles
Definition: mdlj.c:175

Here is the call graph for this function:

Here is the caller graph for this function:

int xyz_in ( FILE *  fp,
double *  rx,
double *  ry,
double *  rz,
double *  vx,
double *  vy,
double *  vz,
int *  N 
)

Definition at line 66 of file mdlj.c.

68  {
69  int i;
70  int has_vel, dum;
71  fscanf(fp,"%i %i\n\n",N,&has_vel);
72  for (i=0;i<(*N);i++) {
73  fscanf(fp,"%i %lf %lf %lf ",&dum,&rx[i],&ry[i],&rz[i]);
74  if (has_vel) { // read velocities
75  fscanf(fp,"%lf %lf %lf",&vx[i],&vy[i],&vz[i]);
76  }
77  }
78  return has_vel;
79 }

Here is the caller graph for this function:

void xyz_out ( FILE *  fp,
double *  rx,
double *  ry,
double *  rz,
double *  vx,
double *  vy,
double *  vz,
int *  ix,
int *  iy,
int *  iz,
double  L,
int  N,
int  z,
int  put_vel,
int  unfold 
)

Definition at line 47 of file mdlj.c.

51  {
52  int i;
53 
54  fprintf(fp,"%i %i\n\n",N,put_vel);
55  for (i=0;i<N;i++) {
56  fprintf(fp,"%i %.8lf %.8lf %.8lf ",z,
57  rx[i]+(unfold?(ix[i]*L):0.0),
58  ry[i]+(unfold?(iy[i]*L):0.0),
59  rz[i]+(unfold?(iz[i]*L):0.0));
60  if (put_vel)
61  fprintf(fp,"%.8lf %.8lf %.8lf",vx[i],vy[i],vz[i]);
62  fprintf(fp,"\n");
63  }
64 }

Here is the caller graph for this function:

Generated on Tue Jan 2 2018 00:06:55 for mdlj by   doxygen 1.8.6