/***********/ /* relax.h */ /***********/ #ifndef _RELAX_H #define _RELAX_H /** universal constants **/ #ifndef PI #define PI 3.14159265358979323846 #endif #define AVO 6.0221367E+23 #define BOLZ 1.380658E-23 /* in Joule */ #define HBAR 1.054589E-34 /* in Joule.s */ /** argon LJ6-12 constants **/ #define SIGMA 3.405E-10 /* in meter */ #define EPSILON (119.8*BOLZ) /* in Joule */ /** reduced units benchmarked to argon LJ6-12 **/ #define UENERGY_IN_JOULE EPSILON #define ULENGTH_IN_METER SIGMA /* from now on are all derived reduced units --> */ #define ULENGTH_IN_ANGSTROM (SIGMA*1E10) #define USTRESS_IN_PASCAL \ (UENERGY_IN_JOULE/ULENGTH_IN_METER/ULENGTH_IN_METER/ULENGTH_IN_METER) #define USTRESS_IN_MPA (USTRESS_IN_PASCAL/1.E6) #define USTRESS_IN_BAR (USTRESS_IN_PASCAL/1.E5) /* Woon's ab initio pair potential */ #define EV_IN_JOULE 1.60217733E-19 #define HARTREE_IN_JOULE 4.3597482E-18 #define BOHR_RADIUS_IN_ANGSTROM 0.529177249 #define WOON_EPSILON (0.4228*HARTREE_IN_JOULE/1000./UENERGY_IN_JOULE) #define WOON_RE (7.1714*BOHR_RADIUS_IN_ANGSTROM/ULENGTH_IN_ANGSTROM) #define WOON_Q 0.6060 #define WOON_KE (0.6627*HARTREE_IN_JOULE/1000./UENERGY_IN_JOULE \ /BOHR_RADIUS_IN_ANGSTROM/BOHR_RADIUS_IN_ANGSTROM \ *ULENGTH_IN_ANGSTROM*ULENGTH_IN_ANGSTROM) #define WOON_K1 sqrt((WOON_Q+1)*WOON_KE/WOON_Q/WOON_EPSILON) #define WOON_K2 sqrt(WOON_Q*WOON_KE/(WOON_Q+1)/WOON_EPSILON) /* D. E. Woon, Chem. Phys. Lett. 204, 29 (1993) */ /* potential cutoff in reduced length unit */ #define CUTOFF (2.5*3.39/ULENGTH_IN_ANGSTROM) #define MAX_STRING_LENGTH 100 /* benchmark state to calculate reference properties */ #define BENCHMARK_NC 4 #define BENCHMARK_NP (4*BENCHMARK_NC*BENCHMARK_NC*BENCHMARK_NC) #define BENCHMARK_LATTICE_CONSTANT (5.3/ULENGTH_IN_ANGSTROM) #define BENCHMARK_EPS (1E-10) double H_benchmark[3][3], stress_benchmark[3][3]; double s_benchmark[3*BENCHMARK_NP], f_benchmark[3*BENCHMARK_NP]; double ep_benchmark, volumep_benchmark; /* provided defects */ const char *provided_defect_names[] = {"perfect_crystal", "vacancy", "interstitial"}; const int provided_defect_excess[] = {0, -1, 1}; int provided_defect_index; /** control variables **/ #define MAX_SHADOW_CAPTURE 256 char defect_name[MAX_STRING_LENGTH], potential_function[MAX_STRING_LENGTH]; int minimize_provided_defect, minimize_user_defined_defect; int LJ6_12, WOON, const_H, const_s, H0_enforced, write_config; int apply_shadow_positioning, use_default_iatom, iatom; double shadow_radius, iatom_s[3]; int save_screen_output; long iseed; /* random number generator seed */ /** common varaibles **/ int np; /* number of particles */ int nc, npfcc; /* reference fcc state */ double H0[3][3], H[3][3], volume, HI[3][3], eta[3][3], M[3][3]; double pote, hamiltonian; /* potential and kinetic energy */ double sout[3][3], stress[3][3]; /* external and internal stress */ double perturbation_amplitude, minimization_tolerance; /** particle variables **/ double *f; /* forces in real and rescaled space */ double *s; /* positions in rescaled space [0,1) */ double *mass; /* particle masses in reduced unit */ /** input/output files **/ char fn_config_read[MAX_STRING_LENGTH], fn_config_write[MAX_STRING_LENGTH]; FILE *config_read, *config_write; char fn_screen_output[MAX_STRING_LENGTH]; FILE *screen_output; /* "tee" screen output to a file */ char fn_xyz_bef[MAX_STRING_LENGTH], fn_xyz_aft[MAX_STRING_LENGTH]; /* Numerical recipes declarations */ #define EPS 1.0e-10 #define GOLD 1.618034 #define GLIMIT 100.0 #define TINY 1.0e-20 #define CGOLD 0.3819660 #define ZEPS 1.0e-10 #define TOL 2.0e-4 #define ITMAX 200 #define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a)) static double maxarg1,maxarg2; #define FMAX(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1) > \ (maxarg2)? (maxarg1) : (maxarg2)) #define MOV3(a,b,c, d,e,f) (a)=(d);(b)=(e);(c)=(f) #define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d) /* intermediate array that f1dim() uses */ static double *f1dim_xt; /* entities that linmin() pass to f1dim() */ static double (*linmin_nrfunc)(double []); static int linmin_ncom; static double *linmin_pcom, *linmin_xicom; /***************************/ /** function declarations **/ /***************************/ double value(double *s); int value_already_called; void force (double *s, double *f); double new_H(double H0[3][3], double eta[3][3], double M[3][3], double H[3][3], double HI[3][3]); void eta_to_M (double eta[3][3], double M[3][3]); double driving_force (double H0[3][3], double eta[3][3], double c); double work_function (double H0[3][3], double eta[3][3], double sout[3][3]); void shadow_positioning (char *fn_xyz); void bubble_sort (int nitems, int *index, double *d); double pair_potential_energy(int np, double H[3][3], double *s); void pair_potential_force (int np, double H[3][3], double *s, double *f, double stress[3][3]); void write_on_config (char filename[]); void read_from_config (char filename[]); void assign_fcc_lattice(int nc, double *s); double frandom(); void mateqv(double A[3][3], double B[3][3]); void matran(double A[3][3], double B[3][3]); double matnorm(double A[3][3]); void matadd(double A[3][3], double B[3][3], double C[3][3]); void matadd_mul(double A[3][3], double CC, double B[3][3], double C[3][3]); void matsub(double A[3][3], double B[3][3], double C[3][3]); double matinv(double A[3][3], double B[3][3]); void matmul(double A[3][3], double B[3][3], double C[3][3]); void polint(double xa[], double ya[], int n, double x, double *y, double *dy); void print(const char *p, ...); /* inaccessible to outside except cg() */ void cg(double p[], int n, double ftol, int *iter, double *fret, double (*func)(double []), void (*dfunc)(double [], double [])); static void linmin (double p[], double xi[], int n, double *fret, double (*func)(double [])); static double f1dim(double x); static void mnbrak (double *ax, double *bx, double *cx, double *fa, double *fb, double *fc, double (*func)(double)); static double brent (double ax, double bx, double cx, double (*f)(double), double tol, double *xmin); static void nrerror(char error_text[]); #endif