| /*************************************************************************/ |
| /* */ |
| /* Copyright (c) 1994 Stanford University */ |
| /* */ |
| /* All rights reserved. */ |
| /* */ |
| /* Permission is given to use, copy, and modify this software for any */ |
| /* non-commercial purpose as long as this copyright notice is not */ |
| /* removed. All other uses, including redistribution in whole or in */ |
| /* part, are forbidden without prior written permission. */ |
| /* */ |
| /* This software is provided with absolutely no warranty and no */ |
| /* support. */ |
| /* */ |
| /*************************************************************************/ |
| |
| EXTERN_ENV |
| |
| #include "math.h" |
| #include "stdio.h" |
| #include "mdvar.h" |
| #include "water.h" |
| #include "cnst.h" |
| #include "fileio.h" |
| #include "parameters.h" |
| #include "mddata.h" |
| #include "split.h" |
| |
| INITIA() |
| { |
| /* this routine initializes the positions of the molecules along |
| a regular cubical lattice, and randomizes the initial velocities of |
| the atoms. The random numbers used in the initialization of velocities |
| are read from the file random.in, which must be in the current working |
| directory in which the program is run */ |
| |
| static double XMIN = 0; |
| static double YMIN = 0; |
| static double ZMIN = 0; |
| FILE *random_numbers; /* points to input file containing |
| pseudo-random numbers for initializing |
| velocities */ |
| double XMAS[4], XS, ZERO, WCOS, WSIN, XT[4], YT[4], Z; |
| double SUX, SUY, SUZ, SUMX, SUMY, SUMZ, FAC; |
| int mol=0; |
| int atom=0; |
| int deriv; |
| double xrand(); |
| |
| random_numbers = fopen("random.in","r"); |
| if (random_numbers == NULL) { |
| fprintf(stderr,"Error in opening file random.in\n"); |
| fflush(stderr); |
| exit(-1); |
| } |
| |
| XMAS[1]=sqrt(OMAS*HMAS); |
| XMAS[0]=HMAS; |
| XMAS[2]=HMAS; |
| |
| /* .....assign positions */ |
| { |
| double NS = pow((double) NMOL, 1.0/3.0) - 0.00001; |
| double XS = BOXL/NS; |
| double ZERO = XS * 0.50; |
| double WCOS = ROH * cos(ANGLE * 0.5); |
| double WSIN = ROH * sin(ANGLE * 0.5); |
| int i,j,k; |
| |
| printf("\nNS = %.16f\n",NS); |
| printf("BOXL = %10f\n",BOXL); |
| printf("CUTOFF = %10f\n",CUTOFF); |
| printf("XS = %10f\n",XS); |
| printf("ZERO = %g\n",ZERO); |
| printf("WCOS = %f\n",WCOS); |
| printf("WSIN = %f\n",WSIN); |
| fflush(stdout); |
| |
| #ifdef RANDOM |
| /* if we want to initialize to a random distribution of displacements |
| for the molecules, rather than a distribution along a regular lattice |
| spaced according to intermolecular distances in water */ |
| srandom(1023); |
| for (i = 0; i < NMOL; i++) { |
| VAR[mol].F[DISP][XDIR][O] = xrand(0, BOXL); |
| VAR[mol].F[DISP][XDIR][H1] = VAR[mol].F[DISP][XDIR][O] + WCOS; |
| VAR[mol].F[DISP][XDIR][H2] = VAR[mol].F[DISP][XDIR][H1]; |
| VAR[mol].F[DISP][YDIR][O] = xrand(0, BOXL); |
| VAR[mol].F[DISP][YDIR][H1] = VAR[mol].F[DISP][YDIR][O] + WSIN; |
| VAR[mol].F[DISP][YDIR][H2] = VAR[mol].F[DISP][YDIR][O] - WSIN; |
| VAR[mol].F[DISP][ZDIR][O] = xrand(0, BOXL); |
| VAR[mol].F[DISP][ZDIR][H1] = VAR[mol].F[DISP][ZDIR][O]; |
| VAR[mol].F[DISP][ZDIR][H2] = VAR[mol].F[DISP][ZDIR][O]; |
| } |
| #else |
| /* not random initial placement, but rather along a regular |
| lattice. This is the default and the prefered initialization |
| since random does not necessarily make sense from the viewpoint |
| of preserving bond distances */ |
| |
| fprintf(six, "***** NEW RUN STARTING FROM REGULAR LATTICE *****\n"); |
| fflush(six); |
| XT[2] = ZERO; |
| mol = 0; |
| for (i=0; i < NS; i+=1) { |
| XT[1]=XT[2]+WCOS; |
| XT[3]=XT[1]; |
| YT[2]=ZERO; |
| for (j=0; j < NS; j+=1) { |
| YT[1]=YT[2]+WSIN; |
| YT[3]=YT[2]-WSIN; |
| Z=ZERO; |
| for (k = 0; k < NS; k++) { |
| for (atom = 0; atom < NATOMS; atom +=1) { |
| VAR[mol].F[DISP][XDIR][atom] = XT[atom+1]; |
| VAR[mol].F[DISP][YDIR][atom] = YT[atom+1]; |
| VAR[mol].F[DISP][ZDIR][atom] = Z; |
| } |
| mol += 1; |
| Z=Z+XS; |
| } |
| YT[2]=YT[2]+XS; |
| } |
| XT[2]=XT[2]+XS; |
| } |
| |
| if (NMOL != mol) { |
| printf("Lattice init error: total mol %d != NMOL %d\n", mol, NMOL); |
| exit(-1); |
| } |
| #endif |
| } |
| |
| /* ASSIGN RANDOM MOMENTA */ |
| fscanf(random_numbers,"%lf",&SUX); |
| |
| SUMX=0.0; |
| SUMY=0.0; |
| SUMZ=0.0; |
| /* read pseudo-random numbers from input file random.in */ |
| for (mol = 0; mol < NMOL; mol++) { |
| for (atom = 0; atom < NATOMS; atom++) { |
| fscanf(random_numbers,"%lf",&VAR[mol].F[VEL][XDIR][atom]); |
| fscanf(random_numbers,"%lf",&VAR[mol].F[VEL][YDIR][atom]); |
| fscanf(random_numbers,"%lf",&VAR[mol].F[VEL][ZDIR][atom]); |
| SUMX = SUMX + VAR[mol].F[VEL][XDIR][atom]; |
| SUMY = SUMY + VAR[mol].F[VEL][YDIR][atom]; |
| SUMZ = SUMZ + VAR[mol].F[VEL][ZDIR][atom]; |
| for (deriv = ACC; deriv < MAXODR; deriv++) { |
| VAR[mol].F[deriv][XDIR][atom] = 0.0; |
| VAR[mol].F[deriv][YDIR][atom] = 0.0; |
| VAR[mol].F[deriv][ZDIR][atom] = 0.0; |
| } |
| } /* atoms */ |
| } /* molecules */ |
| |
| /* find average momenta per atom */ |
| SUMX=SUMX/(NATOMS*NMOL); |
| SUMY=SUMY/(NATOMS*NMOL); |
| SUMZ=SUMZ/(NATOMS*NMOL); |
| |
| /* find normalization factor so that <k.e.>=KT/2 */ |
| SUX=0.0; |
| SUY=0.0; |
| SUZ=0.0; |
| for (mol = 0; mol < NMOL; mol++) { |
| SUX = SUX + (pow( (VAR[mol].F[VEL][XDIR][H1] - SUMX),2.0) |
| +pow( (VAR[mol].F[VEL][XDIR][H2] - SUMX),2.0))/HMAS |
| +pow( (VAR[mol].F[VEL][XDIR][O] - SUMX),2.0)/OMAS; |
| |
| SUY = SUY + (pow( (VAR[mol].F[VEL][YDIR][H1] - SUMY),2.0) |
| +pow( (VAR[mol].F[VEL][YDIR][H2] - SUMY),2.0))/HMAS |
| +pow( (VAR[mol].F[VEL][YDIR][O] - SUMY),2.0)/OMAS; |
| |
| SUZ = SUZ + (pow( (VAR[mol].F[VEL][ZDIR][H1] - SUMZ),2.0) |
| +pow( (VAR[mol].F[VEL][ZDIR][H2] - SUMZ),2.0))/HMAS |
| +pow( (VAR[mol].F[VEL][ZDIR][O] - SUMZ),2.0)/OMAS; |
| } |
| FAC=BOLTZ*TEMP*NATMO/UNITM * pow((UNITT*TSTEP/UNITL),2.0); |
| SUX=sqrt(FAC/SUX); |
| SUY=sqrt(FAC/SUY); |
| SUZ=sqrt(FAC/SUZ); |
| |
| /* normalize individual velocities so that there are no bulk |
| momenta */ |
| XMAS[1]=OMAS; |
| for (mol = 0; mol < NMOL; mol++) { |
| for (atom = 0; atom < NATOMS; atom++) { |
| VAR[mol].F[VEL][XDIR][atom] = ( VAR[mol].F[VEL][XDIR][atom] - |
| SUMX) * SUX/XMAS[atom]; |
| VAR[mol].F[VEL][YDIR][atom] = ( VAR[mol].F[VEL][YDIR][atom] - |
| SUMY) * SUY/XMAS[atom]; |
| VAR[mol].F[VEL][ZDIR][atom] = ( VAR[mol].F[VEL][ZDIR][atom] - |
| SUMZ) * SUZ/XMAS[atom]; |
| } /* for atom */ |
| } /* for mol */ |
| |
| fclose(random_numbers); |
| |
| } /* end of subroutine INITIA */ |
| |
| /* |
| * XRAND: generate floating-point random number. |
| */ |
| |
| double xrand(xl, xh) |
| double xl, xh; /* lower, upper bounds on number */ |
| { |
| long random (); |
| double x; |
| |
| x=(xl + (xh - xl) * ((double) random()) / 2147483647.0); |
| return (x); |
| } |
