cdmc.c

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

/*
************************************************************************
Include Files
************************************************************************
*/
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include "strsub.h"
#include "parse.h"
#include "cdhouse.h"
#include "cdsubs.h"
#include "cdstep.h"


/*
************************************************************************
Global variables
************************************************************************
*/

/*  HALT SIGNAL  */
extern int isig;

/*  GLOBAL VARIABLES  */
extern Simulation_t *s;
extern Particle_t *a;



/*
************************************************************************
Local Function Prototypes
************************************************************************
*/

/*  LOCAL FUNCTION PROTOTYPES  */
double __eharm (Particle_t *, double, double);


/*
************************************************************************
Exported Functions
************************************************************************
*/


/*  perform monte carlo minimization  */

void runmc (Particle_t *a, Simulation_t *s, int nstep, double temp)
	{
	double etmp, eorg, eold, ecur, enew, edel, efac, tfac;
	int j,reject;
	long navgold, nacc=0, istep;
	double esum;
	double *cold = NULL;
	double *cnew = NULL;
	char *tname;
	double eavg=0, evar0=0, evar1=0;
	/*  VARIALBES FOR SAVING COMPARE COORDINATES  */
	FILE *CompareFile;

	/*  INITIALIZE RANDOM JUMPING  */
	int	 rand_half = RAND_MAX >> 1;
	double rsize	  = s->size/rand_half;


	ALLOCATE (cnew, double, 3*a->np)

	/*  Open energy output file  */
	OpenStepOutput ( &(s->EnergyStepOutput) );


	/*  Open Compare file  */
	if (s->UseCompare)
		{

		/*  GET FILE NAME  */
		tname = s->CompareFileName;
		if (*tname=='+')
			tname++;

		/*  OPEN FILE FOR FIRST TIME AS WRITE	*/
		if (s->NewCompareFile)
			{
			s->NewCompareFile = FALSE;
			CompareFile = fopen (tname, "wb");
			}

		/*  OPEN FILE AS APPEND  */
		else
			CompareFile = fopen (tname, "ab");

		if (CompareFile==NULL)
			{
			printf ("ERROR: Cannot open COMPARE file %s.\n",
									s->CompareFileName);
			CleanAfterError();
			}
		}

	/*  USE HARMONIC ENERGY AS "BASE" ENERGY  */
	if (s->UseCompare && s->UseInverseCompare)
		{
		ecur = __eharm (a, s->CompareRadius, s->CompareEnergy);
		}
	/*  USE FULL ENERGY AS "BASE" ENERGY  */
	else
		{
		energy_all (a, s, 0);
		ecur = a->epot;
		}

	/* SAVE OLD PARTICLE POINTER - CHANGE TO NEW PARTICLES  */
	cold	 = a->cur;
	a->cur = cnew;

	navgold = a->navg;
	esum = 0.0;
	if (temp!=0)
		tfac = 1.0/(BK*temp);


	eorg	  = ecur;
	/*  monte carlo minimization	*/
	for (istep=0;istep<nstep && isig!=1;istep++)
		{
		a->step++;
		s->nstep++;

		eold = ecur;
		/*  assign random displacements	*/
		for (j=0;j<3*a->np;j++)
			cnew[j] = cold[j] + rsize * (rand() - rand_half);

		/*  USE HARMONIC ENERGY AS "BASE" ENERGY  */
		if (s->UseCompare && s->UseInverseCompare)
			{
			enew = __eharm (a, s->CompareRadius, s->CompareEnergy);
			}
		/*  USE FULL ENERGY AS "BASE" ENERGY  */
		else
			{
			energy_all (a, s, 0);
			enew = a->epot;
			}
		edel = enew-eold;

		reject = (edel>0);

		if (reject)
			{
			if (temp!=0)
				{
				efac = tfac*edel;
				if (efac<20)
					reject = rand() > RAND_MAX*exp(-efac);
				}
			}


		if (reject)
			{
			s->nrej++;
			ecur = eold;
			}
		else
			{
			for (j=0;j<3*a->np;j++)
				cold[j] = cnew[j];
			ecur = enew;
			s->nacc++;
			nacc++;
			}

		/*  perform running sums  */
		esum += ecur;

		/*  Save energy if saving in effect and step accepted  */
		if (ForceStepOutput ( &(s->EnergyStepOutput), !reject)  )
			fprintf
				(
				s->EnergyStepOutput.File,
				"%li %le\n",
				a->step, s->eunit*ecur/a->np
				);


		/* COMPARE LATTICE  */
		if (!reject && s->UseCompare)
			{
			/*  USE HARMONIC ENERGY AS "COMPARISON" ENERGY  */
			if (!s->UseInverseCompare)
				{
				edel = __eharm (a, s->CompareRadius, s->CompareEnergy) - ecur;
				}
			/*  USE FULL	  ENERGY AS "COMPARISON" ENERGY  */
			else
				{
				energy_all (a, s, 0);
				edel = a->epot - ecur;
				}
			/*  WRITE RESULTS TO FILE	*/
			fwrite (&a->step, sizeof(a->step), 1, CompareFile);
			fwrite (&ecur, 	sizeof(ecur),	  1, CompareFile);
			fwrite (&edel, 	sizeof(edel   ), 1, CompareFile);
			}

		/*  CORRELATION INFO  */
		etmp	 = ecur-eorg;
		evar0 += etmp*etmp;
		evar1 += etmp*(eold-eorg);
		eavg	+= etmp;

		}

	/*  restore particle pointer to original storage location  */
	a->cur = cold;

	if (navgold+a->navg>0)
		a->epavg = (navgold*a->epavg + esum) / (navgold+a->navg);

	FREE (cnew)

	/*  CLOSE FILES  */
	CloseStepOutput ( &(s->EnergyStepOutput) );
	if (s->UseCompare)
		fclose (CompareFile);

	/*  PRINT NUMBER ACCEPTED	*/
	if (istep>0)
		{
		double corr1;

		printf ("Number of Metropolis steps performed %li.\n", istep);
		printf ("Number of Metropolis steps accepted  %li (%6.2f).\n",
					  nacc, (double) nacc / istep);
		eavg	/= istep;
		evar0  = evar0 / istep	- eavg*eavg;
		evar1  = evar1 / istep	- eavg*eavg;
		if (evar0>0)
			corr1 = evar1 / evar0;
		if (evar0<=0.0 || corr1>=1.0)
			printf ("Correlation length: infinite or unknown.\n");
		else
			printf ("Estimated correlation length: %lf\n", 1.0 / (1.0 - corr1) );
		}
	}

void read_compare (char *instr)
	{
	char *tokstr, *fname;
	double *ctemp;

	/*  GET FIRST "OFF" or "ON"  */
	tokstr = strhed (&instr);
	if 	  (!strcmpi(tokstr,"ON"))
		s->UseCompare = TRUE;
	else if (!strcmpi(tokstr,"OFF"))
		{
		s->UseCompare = FALSE;
		return;
		}
	else
		{
		printf ("        ***  ERRROR:  UNKNOWN OPTION");
		return;
		}

	/*  CHECK FOR PRESENCE OF REFERENCE PARTICLES WHEN DOING COMPARE	*/
	if (a->ref == NULL)
			{
		printf ("ERROR (read_compare):  Must have reference particles.\n");
		CleanAfterError();
			}

	/*  CHECK FOR INVERSE OPTION	*/
	tokstr = strhed (&instr);
	if (!strcmpi(tokstr,"INV"))
		{
		s->UseInverseCompare = TRUE;
		tokstr = strhed (&instr);
		}
	else
		s->UseInverseCompare = FALSE;

	/*  READ RADIUS  */
	s->CompareRadius = dblstrf(&tokstr) * 1e-8;

	/*  READ FILE NAME  */
	fname = strhed (&instr);
	strcpy (s->CompareFileName, fname);
	/*  DETERMINE IF FILE OPENED AS WRITE (NOT APPEND) FIRST TIME	*/
	if (*(s->CompareFileName)=='+')
		s->NewCompareFile = FALSE;
	else
		s->NewCompareFile = TRUE;

	/*  GET REFERENCE ENERGY  */
	ctemp = a->cur;
	a->cur = a->ref;
	energy_all (a, s, 0);
	s->CompareEnergy = a->epot;
	a->cur  = ctemp;
	}



/*
************************************************************************
Local Functions
************************************************************************
*/


/*  CALCULATE HARMONIC APPROXIMATION TO ENERGY	*/
double __eharm (Particle_t *a, double CompareRadius, double eref)
	{
	double *ctmp = NULL;
	double  epot, pre_epot, *pre_cur;
	double  *box, box2[3], dmag, t, fac;
	int i,j, n, *surf;
	/*  SET UP BOUNDARY CONDITIONS  */
	surf = a->surf;
	box  = a->bcur;
	for (i=0;i<3;i++)
		box2[i] = box[i]/2;
	/*  CALCULATE MAGNITUDE OF DISPLACEMENT  */
	dmag = 0;
	n = 0;
	for (i=0;i<a->np;i++)
	for (j=0;j<3;j++)
		{
		t = fabs(a->cur[n] - a->ref[n]);
		if (!surf[j])
			if (t>box2[j])
				t = box[j] - t;
		dmag += t*t;
		n++;
		}
	/*  RETURN IF NO DISPLACEMENT RELATIVE TO REFERENCE LATTICE  */
	if (dmag==0.0)
		return (eref);
	/*  ALLOCATE TEMPERORARY COORDINATES  */
	ALLOCATE (ctmp, double, 3*a->np);

	/*  DETERMINE DISPLACEMENT SCALING FACTOR */
	fac  = CompareRadius / sqrt(dmag);
	/*  CALCULATE NEW COORDINATES  */
	n = 0;
	for (i=0;i<a->np;i++)
	for (j=0;j<3;j++)
		{
		t = a->cur[n] - a->ref[n];
		if (!surf[j])
			if (t>box2[j])
				t -= box[j];
			else if (t<-box2[j])
				t += box[j];
		ctmp[n] = a->ref[n] + fac*t;
		n++;
		}
	/*  CALCULATE ENERGY FOR SMALL PERTUBATION  */
	pre_cur	= a->cur;
	pre_epot = a->epot;
	a->cur	= ctmp;
	energy_all (a, s, 0);
	/*  CONVERT TO HARMONIC ENERGY FOR FULL DISPLACEMENT	*/
	epot = eref + (a->epot - eref) / (fac*fac);
	/*  REPLACE ORIGINAL VALUES  */
	a->cur  = pre_cur;
	a->epot = pre_epot;

	/*  FREE STORAGES  */
	FREE (ctmp)

	/*  RETURN ENERGY  */
	return (epot);
	}