cstrain.c

一个很好的分子动力学程序

#include <stdio.h>
#include <strnsub.h>
#include <cor0.h>
#include <math.h>
#include <stdlib.h>

float  *_x;


int _qcmp (const void *i, const void *j) {
   float a = _x[* ( (int *) i) ];
   float b = _x[* ( (int *) j) ];
   if (a<b)
      return(-1);
   if (a>b)
      return( 1);
   return( 0);
}

void _qsorti (int *ix, float *x, int n) {
   _x = x;
   qsort ((void *) ix, n, sizeof (int *), _qcmp);
}


void calc_strain (float radius, stepdata *r, stepdata *s, int selfl)
{
   /* COORDINATES */
   float *rx, *ry, *rz, *cx, *cy, *cz;
/*
   coord  *ref = r->cur;
   coord  *cur = s->cur;
*/
   stepdata *ref = r;
   stepdata *cur = s;
   float cmin[3], cmax[3], cspread[3];
   int n = s->np;
   int ic[3] = {0,1,2};
   int rc[3];
   int *ix, i,j,k, intv, ipart, jpart;
   int bndx, bndy, bndz;
   float t, bx,by,bz,b2x,b2y,b2z;
   float *uxx,*uxy,*uxz,*uyx,*uyy,*uyz,*uzx,*uzy,*uzz,
         *dxx,*dxy,*dxz,*dyy,*dyz,*dzz;
   float txx,txy,txz,tyx,tyy,tyz,tzx,tzy,tzz;
   float qxx,qxy,qxz,qyy,qyz,qzz;
   float q[6];
   float xd,yd,zd, xu,yu,zu, xa,ya,za;
   float rad1,rad2,rad3,radsq, rsqr;
   unsigned *pind, *sind;
   unsigned is,js,ntemp,nalloc;

   /*  CHECK DATA INTEGRITY  */
   if (s->np==0) {
      printf ("ERROR (calc_strain): Number of particles is zero.\n");
      printf ("  Will return to calling program.\n");
      return;
   }
   if (cur->c==NULL || s->t==NULL ||
       cur->c[0]==NULL || cur->c[1]==NULL || cur->c[2]==NULL )  {
      printf ("ERROR (calc_strain) Coordinate arrays for current step not allocated.\n");
      exit(1);
   }
   if (ref->c==NULL || r->t==NULL ||
       ref->c[0]==NULL || ref->c[1]==NULL || ref->c[2]==NULL )  {
      printf ("ERROR (calc_strain) Coordinate arrays for reference step not allocated.\n");
      exit(1);
   }



   /*  ORDER DIRECTIONS BY SPREAD OF cur COORDINATES  */
   for (i=0;i<3;i++)
      cmin[i] = cmax[i] = cur->c[i][0];
   for (i=0;i<3;i++)  {
      for (j=1;j<n;j++)
         if (cur->c[i][j]<cmin[i])
            cmin[i] = cur->c[i][j];
         else if (cur->c[i][j]>cmax[i])
            cmax[i] = cur->c[i][j];
      cspread[i] = cmax[i] - cmin[i];
   }

   /*  SORT (CRUDELY) DIRECTION BY SPREAD  */
   for (i=0;i<2;i++)  {
      for (j=0;j<2;j++)  {
         if (cspread[ic[j]]<cspread[ic[j+1]])  {
            k       = ic[j];
            ic[j]   = ic[j+1];
            ic[j+1] = k;
         }
      }
   }

   /*  SET POINTERS TO DIRECTIONS DERTERMINED BY SORT  */
   cx = cur->c[ic[0]];
   cy = cur->c[ic[1]];
   cz = cur->c[ic[2]];
   rx = (float *) calloc (s->np, sizeof(float) );
   ry = (float *) calloc (s->np, sizeof(float) );
   rz = (float *) calloc (s->np, sizeof(float) );
   if (rz==NULL) {
      printf ("ERROR (calc_strain):  Not enough room to store %i temporary coordinates.\n", s->np);
      exit(1);
   }
   for (i=0;i<s->np;i++) {
      rx[i] = ref->c[ic[0]][i];
      ry[i] = ref->c[ic[1]][i];
      rz[i] = ref->c[ic[2]][i];
   }

   /*  BOX VALUES FOR NEW DIRECTIONS  */
   bndx = s->boundrep[ic[0]];
   bndy = s->boundrep[ic[1]];
   bndz = s->boundrep[ic[2]];
   bx   = cur->box[ic[0]];
   by   = cur->box[ic[1]];
   bz   = cur->box[ic[2]];
   b2x  = cur->box[ic[0]] - radius;
   b2y  = cur->box[ic[1]] - radius;
   b2z  = cur->box[ic[2]] - radius;


   /*  CALCULATE BACKGROUND DISPLACEMENTS  */
   n = r->np;
   for (i=0;i<n;i++)  {
      /*  CALCULATE DISPLACEMENT OF PARTICLES  */
      t = cx[i] - rx[i];
      if (bndx) {
         if (t<= -b2x)
            t +=  bx;
         else if (t>b2x)
            t -=  bx;
      }
      rx[i] = t;

      t = cy[i] - ry[i];
      if (bndy) {
         if (t<= -b2y)
            t +=  by;
         else if (t>b2y)
            t -=  by;
      }
      ry[i] = t;

      t = cz[i] - rz[i];
      if (bndz) {
         if (t<= -b2z)
            t +=  bz;
         else if (t>b2z)
            t -=  bz;
      }
      rz[i] = t;
   }

   /*  CUTOFF RANGES  */
   rad1  = radius*sqrt(1);
   rad2  = radius*sqrt(2);
   rad3  = radius*sqrt(3);
   radsq = radius*radius;

   /*  SORT PARTICLE BY PRIMARY DIRECTION  */
   ix = (int *) calloc (n, sizeof(int) );
   if (ix==NULL) {
      printf ("ERROR (calc_strain):  Not enough room for sorting index\n");
      exit(1);
   }
   for (i=0;i<n;i++)
      ix[i] = i;
   _qsorti (ix, cx, n);


   /*  FIND FIRST ELIGIBLE PARTICLE  */
   if (bndx &&  2*rad1<bx) {
      t    = cx[ix[0]] - rad1 + bx;
      intv = n;
      while (cx[ix[--intv]]>t && intv>0);
      intv++;
      if (intv==n)
         intv=0;
      /* Internal Consistency Check  */
      if (intv<0 || intv>=n) {
         printf ("INTERNAL ERROR (calc_strain): "
                 "intv (%i) out of range [0..%i].\n", intv, n-1);
         printf ("Will reset intv to 0.\n");
         intv = 0;
      }
   }
   else
      intv = 0;


   /*  if selfl create indices into selected particles  */
   if (selfl)
      nalloc = cur->nsel;
   else
      nalloc = cur->np;
   pind = calloc (nalloc,  sizeof(unsigned));
   sind = calloc (cur->np, sizeof(unsigned));
   if (sind==NULL) {
      printf ("ERROR (calc_strain):  Not enough room to store temporary indices.\n");
      exit(1);
   }
   ntemp = 0;
   for (ipart=0;ipart<cur->np;ipart++)
      if (!selfl || cur->sel[ipart]) {
         if (ntemp<nalloc)
            pind[ntemp] = ipart;
         sind[ipart] = ntemp;
         ntemp++;
      }
   if (ntemp>nalloc) {
      printf ("*** INTERNAL ERROR (calc_strain): ntemp %i exceeds expected %i\n",
              ntemp, nalloc);
      exit(1);
   }

   /*  CREATE SPACE FOR PARTICLE STATISTICS  */
   uxx = (float *) calloc (nalloc, sizeof(float) );
   uxy = (float *) calloc (nalloc, sizeof(float) );
   uxz = (float *) calloc (nalloc, sizeof(float) );
   uyx = (float *) calloc (nalloc, sizeof(float) );
   uyy = (float *) calloc (nalloc, sizeof(float) );
   uyz = (float *) calloc (nalloc, sizeof(float) );
   uzx = (float *) calloc (nalloc, sizeof(float) );
   uzy = (float *) calloc (nalloc, sizeof(float) );
   uzz = (float *) calloc (nalloc, sizeof(float) );
   dxx = (float *) calloc (nalloc, sizeof(float) );
   dxy = (float *) calloc (nalloc, sizeof(float) );
   dxz = (float *) calloc (nalloc, sizeof(float) );
   dyy = (float *) calloc (nalloc, sizeof(float) );
   dyz = (float *) calloc (nalloc, sizeof(float) );
   dzz = (float *) calloc (nalloc, sizeof(float) );
   if (dzz==NULL) {
      printf ("\nERROR (calc_strain):  Not enough memory to perform strain calculation.\n");
      printf ("Try again with fewer particles.\n");
      printf ("Current number of particles is %i.\n\n", n);
      exit(1);
   }
   /*  START NEIGHBOR SEARCH  */
   for (ipart=0;ipart<n;ipart++) {
      i  = ix[ipart];
      is = sind[i];
      /*  SEARCH NEIGHBORS  */
      jpart = intv;
      while (jpart!=ipart) {
         j  = ix[jpart];
         js = sind[j];
         xd = cx[i]-cx[j];
         if (bndx)
            if (xd>b2x)
               xd -= bx;
            else if (xd<-b2x)
               xd += bx;
         if ((xa=fabs(xd))>=rad1)
            intv = jpart;
         else if (!selfl || s->sel[i] || s->sel[j]) {
            yd = cy[i]-cy[j];
            if (bndy)
               if (yd>b2y)
                  yd -= by;
               else if (yd<-b2y)
                  yd += by;
            if ((ya=fabs(yd))<rad1)
            if (xa+ya<rad2) {
               zd = cz[i]-cz[j];
                  if (bndz) {
                  if (zd>b2z)
                     zd -= bz;
                  else if (zd<-b2z)
                     zd += bz;
               }
               if ((za=fabs(zd))<rad1)
               if (xa+ya+za<rad3) {
                 rsqr =  xd*xd + yd*yd + zd*zd;
                 if (rsqr<radsq)  {
                    xu  = rx[i] - rx[j];
                    yu  = ry[i] - ry[j];
                    zu  = rz[i] - rz[j];
                    if (!selfl || s->sel[i]) {
							  uxx[is] += xu*xd;
                       uxy[is] += xu*yd;
                       uxz[is] += xu*zd;
                       uyx[is] += yu*xd;
                       uyy[is] += yu*yd;
                       uyz[is] += yu*zd;
                       uzx[is] += zu*xd;
                       uzy[is] += zu*yd;
                       uzz[is] += zu*zd;
                       dxx[is] += xd*xd;
                       dxy[is] += xd*yd;
                       dxz[is] += xd*zd;
                       dyy[is] += yd*yd;
                       dyz[is] += yd*zd;
                       dzz[is] += zd*zd;
                    }
                    if (!selfl || s->sel[j]) {
                       uxx[js] += xu*xd;
                       uxy[js] += xu*yd;
                       uxz[js] += xu*zd;
                       uyx[js] += yu*xd;
                       uyy[js] += yu*yd;
                       uyz[js] += yu*zd;
                       uzx[js] += zu*xd;
                       uzy[js] += zu*yd;
                       uzz[js] += zu*zd;
                       dxx[js] += xd*xd;
                       dxy[js] += xd*yd;
                       dxz[js] += xd*zd;
                       dyy[js] += yd*yd;
                       dyz[js] += yd*zd;
                       dzz[js] += zd*zd;
                    }
                 }
               }
            }
         }
         jpart++;
         if (jpart>=n)
            jpart -= n;
      }
   }
   /*  rc ARRAY MAPS COORDINATES BACK INTO ORIGINAL XYZ  */
   for (i=0;i<3;i++)
      if (ic[i]<0 || ic[i]>2)
      {
         printf ("INTERNAL ERROR (calc_strain): direction indices out of range.\n");
         exit(1);
      }
      else
         rc[ic[i]] = i;
   /*  CALCULATE LOCAL STRAIN MATRIX  */
   printf ("BOX %f %f %f\n", cur->box[0], cur->box[1], cur->box[2]);
   printf ("STRAIN %i\n", nalloc);
   for (i=0;i<nalloc;i++)
      if (!selfl || s->sel[pind[i]])  {
         qxx = dxx[i];
         qxy = dxy[i];
         qxz = dxz[i];
         qyy = dyy[i];
         qyz = dyz[i];
         qzz = dzz[i];
         inverse (&qxx,&qxy,&qxz,&qyy,&qyz,&qzz);
         txx = uxx[i] * qxx  +  uxy[i] * qxy  +  uxz[i] * qxz;
         txy = uxx[i] * qxy  +  uxy[i] * qyy  +  uxz[i] * qyz;
         txz = uxx[i] * qxz  +  uxy[i] * qyz  +  uxz[i] * qzz;
         tyx = uyx[i] * qxx  +  uyy[i] * qxy  +  uyz[i] * qxz;
         tyy = uyx[i] * qxy  +  uyy[i] * qyy  +  uyz[i] * qyz;
         tyz = uyx[i] * qxz  +  uyy[i] * qyz  +  uyz[i] * qzz;
         tzx = uzx[i] * qxx  +  uzy[i] * qxy  +  uzz[i] * qxz;
         tzy = uzx[i] * qxy  +  uzy[i] * qyy  +  uzz[i] * qyz;
         tzz = uzx[i] * qxz  +  uzy[i] * qyz  +  uzz[i] * qzz;
         /*  convert from displacment matrix to transformation matrix      */
         /*    displacement   matrix: operates on positions to yield disp. */
         /*    transformation matrix: operates on positions to yield pos.  */
         txx = txx + 1;
         tyy = tyy + 1;
         tzz = tzz + 1;
         /*  take square of t  */
         qxx = txx * txx  +  tyx * tyx  +  tzx * tzx;
         qxy = txx * txy  +  tyx * tyy  +  tzx * tzy;
         qxz = txx * txz  +  tyx * tyz  +  tzx * tzz;
         qyy = txy * txy  +  tyy * tyy  +  tzy * tzy;
         qyz = txy * txz  +  tyy * tyz  +  tzy * tzz;
         qzz = txz * txz  +  tyz * tyz  +  tzz * tzz;
         /*  find pure finite strain  */
         pure (&qxx,&qxy,&qxz,&qyy,&qyz,&qzz);
         /*  convert back to displacement matrix (strain)  */
         qxx = qxx - 1;
         qyy = qyy - 1;
         qzz = qzz - 1;
         /*  map coordinates back into original XYZ  */
         q[0] = qxx; q[1] = qyy; q[2] = qzz;
         q[3] = qyz; q[4] = qxz; q[5] = qxy;
			qxx = q[rc[0]];
         qyy = q[rc[1]];
         qzz = q[rc[2]];
         qyz = 2*q[rc[0]+3];
         qxz = 2*q[rc[1]+3];
         qxy = 2*q[rc[2]+3];
         /*  write type and coordinate   */
         ipart = pind[i];
         printf 
				(
				"%i %8.4f %8.4f %8.4f", 
				s->t[ipart]+1,
				1e8*cur->c[X][ipart], 
				1e8*cur->c[Y][ipart], 
				1e8*cur->c[Z][ipart]
				);
         printf (" %+f %+f %+f %+f %+f %+f\n",
                   qxx, qyy, qzz, qyz, qxz, qxy);
		}
   /*  RELEASE MEMORY  */
   free (dzz);
   free (dyz);
   free (dyy);
   free (dxz);
   free (dxy);
	free (dxx);
   free (uzz);
   free (uzy);
   free (uzx);
   free (uyz);
	free (uyy);
   free (uyx);
   free (uxz);
   free (uxy);
   free (uxx);
   free (pind);
   free (sind);
	free (ix );
   free (rz );
   free (ry );
   free (rx );
}