Generated from cdfastrlx.c with ROBODoc v3.2.2 on Fri Sep 14 14:23:14 2001
TABLE OF CONTENTS
- QMD-MGDFT/cdfastrlx.c
NAME
Ab initio real space code with multigrid acceleration
Quantum molecular dynamics package.
Version: 2.1.5
COPYRIGHT
Copyright (C) 1995 Emil Briggs
Copyright (C) 1998 Emil Briggs, Charles Brabec, Mark Wensell,
Dan Sullivan, Chris Rapcewicz, Jerzy Bernholc
Copyright (C) 2001 Emil Briggs, Wenchang Lu,
Marco Buongiorno Nardelli,Charles Brabec,
Mark Wensell,Dan Sullivan, Chris Rapcewicz,
Jerzy Bernholc
FUNCTION
void cdfastrlx(STATE *states, REAL *vxc, REAL *vh, REAL *vnuc,
REAL *rho, REAL *rhocore, REAL *rhoc)
drive routine for constrained fast relax.
INPUTS
states: all wave functions (see md.h)
vxc: exchange-correlation potential
vh: Hartree potential
vnuc: pseudopotential
rho: total charge density
rhocore: charge density of core electrons, only useful when we
include non-linear core correction for pseudopotential.
rhoc: compensating charge density
OUTPUT
all the above inputs are updated
PARENTS
md.c
CHILDREN
quench.c to_crystal.c to_cartesian.c init_nuc.c get_nlop.c scf.c
get_te.c force.c md_fastrelax.c
SOURCE
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include "md.h"
/* Local function prototypes */
void md_fastrelax(void);
void movie(FILE *);
void cdfastrlx(STATE *states, REAL *vxc, REAL *vh, REAL *vnuc,
REAL *rho, REAL *rhocore, REAL *rhoc)
{
int steps, mdsteps, it, ion;
int CONVERGENCE = FALSE;
int CONV_FORCE = FALSE;
int i;
int temp_num;
int diagcount = 0;
REAL length = 0.0;
REAL projection = 0.0;
REAL vector[3];
ION *iptr;
FILE *mfp=NULL;
if (ct.override_atoms==1) {
quench(states, vxc, vh, vnuc, rho, rhocore, rhoc);
}
/* open movie file */
if(pct.thispe == 0) {
mfp = fopen("traj.rmv","w");
if(setvbuf(mfp,(char *)NULL,_IOFBF,4096*16) != 0)
printf("\n Warning: cant allocate movie io buffer size\n");
/* output initial frame */
movie(mfp);
}
/* normalize direction vector for constrained dynamics */
for (i = 0; i < 3; i++){
length += ct.cd_vector[i]*ct.cd_vector[i];
} /* i */
length = sqrt(length);
if (length <= 0.1e-12) error_handler("cdfastrlx", "null vector");
for (i = 0; i < 3; i++){
vector[i] = ct.cd_vector[i]/length;
}
CONV_FORCE = TRUE;
for(ion = 0;ion < ct.num_ions;ion++) {
if(ct.ions[ion].movable){
CONV_FORCE = (CONV_FORCE &&
((fabs(ct.ions[ion].force[ct.fpt[0]][0]) < ct.thr_frc) &&
(fabs(ct.ions[ion].force[ct.fpt[0]][1]) < ct.thr_frc) &&
(fabs(ct.ions[ion].force[ct.fpt[0]][2]) < ct.thr_frc)));
} /* end if */
} /* ion */
/* relax the system along the non-constrained directions */
for(mdsteps = 0;mdsteps < ct.num_steps;mdsteps++) {
/* enforce periodic boundary conditions on the ions */
for(it = 0;it < ct.num_ions;it++){
iptr = &ct.ions[it];
to_crystal(iptr->xtal,iptr->crds);
if(iptr->xtal[0] > ONE) iptr->xtal[0] -= ONE;
if(iptr->xtal[0] < ZERO) iptr->xtal[0] += ONE;
if(iptr->xtal[1] > ONE) iptr->xtal[1] -= ONE;
if(iptr->xtal[1] < 0.0) iptr->xtal[1] += ONE;
if(iptr->xtal[2] > ONE) iptr->xtal[2] -= ONE;
if(iptr->xtal[2] < 0.0) iptr->xtal[2] += ONE;
to_cartesian(iptr->xtal,iptr->crds);
} /* end for */
if (CONV_FORCE) {
for(ion = 0;ion < ct.num_ions;ion++) {
if(ct.ions[ion].movable==2){
if(pct.thispe == 0) {
printf("\n\n FORCES CONVERGED");
printf("\n updating constrained ion position");
} /* end print */
ct.ions[ion].crds[0] += ct.cd_velocity * vector[0];
ct.ions[ion].crds[1] += ct.cd_velocity * vector[1];
ct.ions[ion].crds[2] += ct.cd_velocity * vector[2];
} /* end if */
/* check if the projection of the change of the constrained
ion's coordinate along the normal direction is longer than
the length of the vector input, i.e. has the ion moved too
far? */
projection = 0.0;
length = 0.0;
for (i = 0; i < 3; i++){
projection += (ct.ions[ion].crds[i]-ct.ions[ion].icrds[i])*vector[i];
length += ct.cd_vector[i]*ct.cd_vector[i];
} /* i */
length = sqrt(length);
if (projection > length) {
if(pct.thispe == 0) {
printf("\n\n constrained atom moved beyond desired endpoint");
printf("\n STOPPING");
return;
}
} /* end if */
} /* ion */
/* Update items that change to incorporate the change of the ion*/
init_nuc(vnuc, rhoc, rhocore);
get_nlop();
quench(states, vxc, vh, vnuc, rho, rhocore, rhoc);
} /* end if CONV_FORCE */
CONVERGENCE = FALSE;
for(steps = 0;steps < ct.scfpermd;steps++) {
/* Perform a single self-consistent step */
if (! CONVERGENCE) {
scf(states, vxc, vh, vnuc, rho, rhocore, rhoc, &CONVERGENCE);
/* Do diagonalizations if requested */
if(ct.diag) {
if(diagcount > ct.diag) {
subdiag(states, vh, vnuc, vxc);
diagcount = 0;
} /* end if */
diagcount++;
} /* end if */
/* Get the total energy */
get_te(rho, rhocore, rhoc, vh, vxc, states);
} /* end if "convergence" */
} /* end for steps */
/* Compute the forces */
force(rho, rhoc, vh, vxc, states);
/* Project out the component of the force along the constraint direction
for the constrained atom */
for(ion = 0;ion < ct.num_ions;ion++) {
if(ct.ions[ion].movable==2){
length = 0.0;
for (i = 0; i < 3; i++){
length += ct.ions[ion].force[ct.fpt[0]][i] * vector[i];
}
ct.ions[ion].force[ct.fpt[0]][0] -= length * vector[0];
ct.ions[ion].force[ct.fpt[0]][1] -= length * vector[1];
ct.ions[ion].force[ct.fpt[0]][2] -= length * vector[2];
} /* end if */
} /* ion */
/* Output the forces */
if(pct.thispe == 0)write_force();
CONV_FORCE = TRUE;
for(ion = 0;ion < ct.num_ions;ion++) {
if(ct.ions[ion].movable){
CONV_FORCE = (CONV_FORCE &&
( (fabs(ct.ions[ion].force[ct.fpt[0]][0]) < ct.thr_frc) &&
(fabs(ct.ions[ion].force[ct.fpt[0]][1]) < ct.thr_frc) &&
(fabs(ct.ions[ion].force[ct.fpt[0]][2]) < ct.thr_frc)));
} /* end if "movable ions" */
} /* end for ion */
if (! CONV_FORCE) {
/* Move the ions */
md_fastrelax();
/* Update items that change when the ionic coordinates change */
init_nuc(vnuc, rhoc, rhocore);
get_nlop();
/* write out frame */
if(pct.thispe == 0) {
movie(mfp);
}
}
} /* end for */
/* close moviefile */
if(pct.thispe == 0) {
fclose(mfp);
}
} /* end fastrlx */