Program Code: Source Files for PR Bulk Simulations
Program Code: Source Files for PR Bulk Simulations (old)
Program Code: Source Files for NVT Bulk Simulations
Program Code: Source Files for Surface Simulations
The MD program ( for surface premelting simulations) accepts an input file ``initdata.h'' where the
size of the computational cell and all the relevant parameters
(which are not changed frequently) are assigned.
//************** The input control file: ``initdata.h''************//
#ifndef INITDATA_H
#define INITDATA_H
#define surface 0 // 0- [001], 1-
[011], 2-[111]
# if surface != 1
#define size_x 10 // number of x-size atoms
#define size_y 10 // number of y-size atoms
#define Atoms 2700 // Number of atoms =size_x*size_y*size_z
#define FROZEN_ATOMS 300 // NUMBER_OF_FROZEN_LAYERS*size_x*size_y
# else //[011]
#define size_x 7 // number of size atoms
#define size_y 7 // number of size atoms
#define Atoms 2646 // Number of atoms =2*size_x*size_y*size_z
#define FROZEN_ATOMS 294 // 2*NUMBER_OF_FROZEN_LAYERS*size_x*size_y
#endif
#define size_z 27 // number of size atoms
#define INIT_LAYERS size_z
#define NUMBER_OF_FROZEN_LAYERS 3
#define NUMBER_OF_LAYERS INIT_LAYERS-NUMBER_OF_FROZEN_LAYERS+1
#define a00 3.0399 /* size of initail cell at T=0 K*/
#define alpha_t 0.0000086 /*Thermal Expansion Coefficient*/
#define dt 0.002 /*Time step*/
#define N_neigh_step 100 /*Number of neighbour steps*/
#define step_EPT_meas 30000
#define step_backup 200000
//Pressure & Temperature Control
#define press 0.0
// Density Profiles
#define DELTA 20
#define RESOLUTION DELTA*(INIT_LAYERS+5)
//Chen_density
# define Z_RES 500
// Pair_Correlation Data
# define PAIR_STEP 200
// Order Parameter
# define NUMBER_OF_DIRECTIONS 3
// Diffusion Date
# define TMAX 3000 /* number of measurements*/
# define T0max 300 /* max number of t0 time origins*/
# define NUMBER_OF_CALLS 10
/* time interval between set of a new t0 time origin*/
#endif
An example of the input file ``aaa_info.txt''
where the parameters which are most frequently changed
(type of surface, temperature, number of measurements, etc)
are assigned. This file is generated automatically at the first run.
***************''aaa_info.txt'' *****************
2
1800
1000
100
100
30000
100
100.0
6.000000
2700
100
3
24
[001]
*************** Main Parameters *****************
Regime: 0- start, 1 - equilibration, 2- measurements
Temperature
Number of steps in order to approach to the equilibrium
Number of steps between two succesive measurements
Number of measurements
Initial step
Number of the aviz pictures
Percent of the accomplished job
NVE/NVT parameter
Number of the atoms
Number of the atoms in a layer
Number of frozen layers
Number of free layers
The surface type
*************** `` func_cp.h '' *****************
The `` func_cp.h '' file describes all the functions are used by the main ``md.c'' program.
# ifndef NPTFUNC_H
# define NPTFUNC_H
void add_info();
double scale_lattice();
void aviz_layer(char *fn);
void check_initial_data();
void layer_colors(void);
void comments(char *fn);
long get_inform(char *fn);
void save_inform(char *fn);
void backup(char *fn, int mode);
void save_configuration(int mode);
void equilibration();
void measurements(void);
void measurements_initialize_n_backup(int mode);
void measurements_make_n_output(int mode);
// Trajectories
void init_trajectories(int mode);
void trajectories();
// Layer distribution
void energy_temp_pressure();
void atom_distribution_in_layers (int mode);
void atom_distribution_in_layers_rw (char *fn, int mode);
//Diffusion
void diffusion_rw(char *fn,int mode);
void find_diffusion_times(void);
void diffusion(int mode);
void pbc_rw(char *fn,int mode);
// Pair Correlation Functions
void pair_correlation(int mode);
void pair_correlation_rw(char *fn, int mode );
// Inter Layer Distances
void inter_layer_dist(int mode);
void inter_layer_dist_rw(char *fn,int mode);
// Order Profile Functions
void order_parameters(int mode);
void order_parameters_rw(char *fn, int mode );
// Temperature Profile Functions
void temperature_profile (int mode);
void temperature_profile_rw(char *fn, int mode);
// Density Profile Functions
void density_profile(int mode);
void density_profile_rw(char *fn, int mode);
void chen_density_profile(int mode);
void chen_density_profile_rw(char *fn, int mode);
//Initalization
void init_r_init(void);
void init(void);
void aviz(char *fn);
//Setup
void bcc_011(void) ;
void bcc_001(void);
void velocity_init(void);
double lattice_parameter(double x);
void setup_111();
void create_001();
//Input Output
void read_coord_from_file(char *);
void read_vel_from_file(char*);
void save_to_file(long,int);
// Basis
void build_neigh_list(void);
void predictor(void);
void corrector(void);
void compute_f(void );
void eval_PTO(void);
// Density part of potential
double density(double pr_dis);
double density_dot(double pr_dis);
double density_dot_dot(double pr_dis);
double density_cut(double r);
double density_dot_cut(double r);
// Phi potential
double phi_pot(double pr_dis);
double phi_pot_dot(double pr_dis);
double phi_pot_dot_dot(double pr_dis);
double phi_pot_cut(double r);
double phi_pot_dot_cut(double r);
double phi_pot_dot_ren(double x);
double phi_pot_dot_dot_ren(double x);
//F_potential
double F_pot(double);
double F_pot_dot(double);
double F_pot_dot_dot(double x);
double max3_plus(double R1,double R2,double R3,double M);
#endif
*************** ``potential.h'' *****************
// All the parameters of the modified FS potential for vanadium are collected in this file
# ifndef POTENTIAL_H
# define POTENTIAL_H
#define cutoff 3.8 /*Cutoff of potential*/
#define neicutoff 14.440000 /*Square of the nearest neighbour distance */
# define FS_d 3.692767
# define r_cut 3.8
# define c1 -0.8816318
# define c2 1.4907756
# define c3 -0.3976370
# define A -2.010637
# define FND 2.63185120210091
# define B 23.0
# define alp 0.5
# define b0 2.632
#endif
*************** ``constants.h'' *****************
The various constants (
,
) are defined here
# ifndef CONSTANT_H
# define CONSTANT_H
#define pi 3.1415926535897932384626433832795
#define kb 0.0000862 /* Boltzman constant */
#endif
*************** ``md.c'' *****************
The main program ``md.c'' which calls other functions.
# include stdio.h
# include stdlib.h
# include math.h
# include "func_cp.h"
# include "potential.h"
# include "initdata.h"
# include "constants.h"
int Color[Atoms ];
double r[Atoms][3];
double v[Atoms][3];
double f[Atoms][3];
int nei[Atoms][Atoms];
// Predictor-Corrector
double f1[Atoms][3];
double f2[Atoms][3];
double r1[Atoms][3];
double v1[Atoms][3];
//General properties
double a0,Q,temp,ksi=0;
double temp_current,press_current;
double pot_E,kin_E,virsum;
double lx,ly,lz,Width;
double pten[3][3];
int regime;
//Layer structure
int atoms_distribution[NUMBER_OF_LAYERS];
long counter_of_atoms_distribution;
long step0,step,num_eq_steps,num_of_meas;
int interval,Limit_aviz_pict,number_aviz_pict=0,step_Aviz;
int layers[NUMBER_OF_LAYERS][Atoms ];
// Pair Correlation Function
double Profile[NUMBER_OF_LAYERS][PAIR_STEP],width_rdf,
normFac,half_box;
long pair_corr_count;
// Order parameter profile
double S[NUMBER_OF_LAYERS][NUMBER_OF_DIRECTIONS][2];
// Chen Density profile
int z_profile[RESOLUTION];
long count_density_profile ;
double dens_ch[Z_RES],temp_ch[Z_RES], chen_counter;
//Diffusion
long diff_counter,tmax,t0max;
int ntime[TMAX],time0[T0max],t0=0;
double r0[Atoms][3][T0max];
double r2f[TMAX][3][NUMBER_OF_LAYERS];
int pbc[Atoms][2];
//Temperature profile
double dis_bet_lay[NUMBER_OF_LAYERS], counter_dist;
double temp_profile[NUMBER_OF_LAYERS], counter_temp;
int main()
{
int i,j,dd;
FILE* outd;
char fn[50];
check_initial_data();
get_inform("aaa_info.txt");
// Take thermal expansion into account
a0=lattice_parameter(temp);
if(surface==0)//==== Structure b.c.c. [001] ==
{
bcc_001() ;
Width=a0/2;
lx=a0*size_x; // a0 Angstrems lattice constants
ly=a0*size_y; // lx, ly length of the box Angstrems lattice constants
} //=================
if(surface==1)
{
Width=a0/sqrt(2.); // Distance between layers in ordered crystal
lx=sqrt(2.)*a0*size_x; // a0 Angstrems lattice constants
ly=a0*size_y; // lx, ly length of the box Angstrems lattice constants
bcc_011() ;
}
//================
if(surface==2)
{
Width=a0/(2*sqrt(3.));
lx=sqrt(2.)*a0*size_x; // a0 Angstrems lattice constants
ly=sqrt(3./2.)*a0*size_y; // lx, ly length of the box Angstrems lattice constants
read_coord_from_file("r00.log");
scale_lattice();
}
//=========================//
lz=Width*(INIT_LAYERS-NUMBER_OF_FROZEN_LAYERS);
layer_colors();
aviz("van_0.xyz");
add_info();
//====================//
step=step0;
init();
if(regime){ read_coord_from_file("r.log"); read_vel_from_file("v.log");}
if(regime 2)
{
if(regime==0)
{
velocity_init();
aviz("van_0.xyz");
regime=1;
}
else fprintf(stderr,"\n Continue approaching to the equilibrium");
//Continue simulation without the structure measurements
equilibration();
step0=step;
regime=2;
}
if(regime==2)
{
fprintf(stderr,"\n Start measurements");
measurements();
}
fprintf(stderr,"\n\t ********* The Happy End ***************\n");
return(0);
}
*************** ``initialization.c'' *****************
Set to zero all the parameters of the system
# include stdio.h
# include stdlib.h
# include "initdata.h"
void init(void)
{
extern double f1[Atoms][3];
extern double f2[Atoms][3];
extern double r1[Atoms][3];
extern double v1[Atoms][3];
extern int z_profile[RESOLUTION];
int i,j;
for(i=0;iAtoms;i++)
for(j=0;j3;j++)
r1[i][j]=v1[i][j]=f1[i][j]=f2[i][j]=0.;
for(i=0;iRESOLUTION;i++)
z_profile[i]=0;
}
///////////////////////////////////////////////////////////////
void pbc_rw(char *fn,int mode)
{
extern int pbc[Atoms][2];
extern double r[Atoms][3];
FILE *fd;
int i;
if(mode==0)
{
///////////////////////////////////////////
if((fd=fopen(fn,"r"))!=NULL)
{
for (i=FROZEN_ATOMS;iAtoms;i++)
fscanf(fd," \n %d %d ",&pbc[i][0],&pbc[i][1]);
fclose(fd);
}
else
{
fprintf(stderr,"\n ! pbc initialization (:) ");
system("rm -rf diffusion.back atom*.txt ");
if((fd=fopen("r0.log","w"))!=NULL)
{
for (i=0;iAtoms;i++)
{
fprintf(fd,"\n %lf %lf %lf", r[i][0],r[i][1],r[i][2]);
pbc[i][0]=pbc[i][1]=0;//Check & keep p.b.c. for each atom
}
fclose(fd);
}
}
//////////////////////////////////////////////
}
if(mode==1)
{
///////////////////////////////////////////
if((fd=fopen(fn,"w"))!=NULL)
{
for (i=FROZEN_ATOMS;iAtoms;i++)
fprintf(fd," \n %d %d ",pbc[i][0],pbc[i][1]);
fclose(fd);
}
}
}
*************** ``equilibration.c'' *****************
Equilibration of the system
# include stdio.h
# include stdlib.h
# include math.h
# include "func_cp.h"
# include "potential.h"
# include "initdata.h"
# include "constants.h"
void equilibration()
{
extern int number_aviz_pict,Limit_aviz_pict,step_Aviz;
extern long step,num_eq_steps;
extern double temp_current,press_current,kin_E,pot_E;
int i,j,eq_steps;
char fn[50];
eq_steps=(int)1.0*num_eq_steps/N_neigh_step;
build_neigh_list();
compute_f();
//Launching calculations
for(i=0;ieq_steps;i++){
for(j=0;jN_neigh_step;j++)
{
step++;
predictor();
compute_f();
eval_PTO();
corrector();
}
build_neigh_list();
fprintf(stderr,"\n step=%ld temp=%g pressure=%g
total_E=%g ",step,temp_current,press_current,(kin_E+pot_E));
if(step % step_EPT_meas==0) energy_temp_pressure();
if(step % step_Aviz==0)
{
if(number_aviz_pict++ Limit_aviz_pict)
{ sprintf(fn,"van_%ld.xyz",step);
aviz(fn);
}
}
if(step % step_backup==0)backup("aaa_info.txt",1);
} // The Main Cycle: over i
backup("aaa_info.txt",1);
}
*************** ``measurements.c'' *****************
Measurement of the various physical properties
# include stdio.h
# include stdlib.h
# include math.h
# include "func_cp.h"
# include "potential.h"
# include "initdata.h"
# include "constants.h"
void measurements()
{
extern long step,num_of_meas,step_Aviz;
extern int interval,Limit_aviz_pict,number_aviz_pict;
extern double temp_current,press_current;
extern double pot_E,kin_E;
long number_of_steps;
int i, j, write_intermediate_results;
char fn[50];
write_intermediate_results=(int)1.0*num_of_meas*interval/4.0;
number_of_steps=(long)1.*num_of_meas*interval/N_neigh_step;
measurements_initialize_n_backup(0); // Initialize or read some functions:
// Starting calculations
build_neigh_list();
compute_f();
for(i=0;inumber_of_steps;i++){
//====================================//
for(j=0;jN_neigh_step;j++)
{
step++;
predictor();
compute_f();
eval_PTO();
corrector();
}
build_neigh_list();
fprintf(stderr,"\n step=%ld temp=%g pressure=%g
total_E=%g ",step,temp_current,press_current,(kin_E+pot_E));
if(step % interval==0)
{
if(step % write_intermediate_results) measurements_make_n_output(1);
else measurements_make_n_output(2) ;
}
if(step % step_Aviz==0)
{
if(number_aviz_pict++ Limit_aviz_pict)
{
sprintf(fn,"van_%ld.xyz",step);
aviz(fn); aviz_layer(fn);
}
}
if(step % step_backup==0) measurements_initialize_n_backup(1);
//=========================//
} // Main cycle
// Final measurements:
measurements_make_n_output(2) ;
// Final backup:
measurements_initialize_n_backup(1);
}
////////////////////////
void measurements_initialize_n_backup(int mode)
{
backup("aaa_info.txt",mode);
density_profile_rw("dens_prof.back",mode);
chen_density_profile_rw("chen_dens_prof.back",mode);
atom_distribution_in_layers_rw ("at_dist_prof.back",mode);
temperature_profile_rw("temp_prof.back",mode);
order_parameters_rw("order_prof.back",mode);
inter_layer_dist_rw("inter_layer_dist.back",mode);
pair_correlation_rw("pair_func.back",mode);
diffusion_rw("diffusion.back",mode);
pbc_rw("pbc.back",mode);
}
///////////////////////////////
void measurements_make_n_output(int mode)
{
energy_temp_pressure();
chen_density_profile(mode);
atom_distribution_in_layers(mode);
temperature_profile(mode);
order_parameters(mode);
inter_layer_dist(mode);
pair_correlation(mode);
diffusion(mode);
}
*************** ``bcc_setup.c'' *****************
How to construct the various low index faces of a bcc crystal
and initialize the velocities
# include stdio.h
# include math.h
#include stdlib.h
#include time.h
# include "initdata.h"
# include "constants.h"
# include "potential.h"
void velocity_init(void){
extern double v[Atoms][3],temp;
int j,k;
double sum[3];// Initial Temperature
double x1,x2,x3,y1,y2,y3,Sf;
for(j=0;jFROZEN_ATOMS;j++)
{ v[j][0]=0.; v[j][1]=0.; v[j][2]=0.; }
for (j=0;j3;j++)
sum[j]=0.;
srand( (unsigned)time( NULL ) );
/* Seed the random-number generator with current
time so that * the numbers will be different every time we run */
/* Generate gaussian distribution */
for (j=FROZEN_ATOMS;jAtoms;j++)
{
x1=((double)rand())/RAND_MAX;
x2=((double)rand())/RAND_MAX;
x3=((double)rand())/RAND_MAX;
/*The Box Muller method */
y1=sqrt(-log(x1))*cos(2*pi*x2);
y2=sqrt(-log(x2))*cos(2*pi*x1);
y3=sqrt(-log(x3))*cos(2*pi*x1);
v[j][0]=y1;
v[j][1]=y2;
v[j][2]=y3;
for (k=0;k3;k++)sum[k]+=v[j][k];
}
/* Center of mass velocity */
for (k=0;k3;k++)
sum[k]/=(Atoms-FROZEN_ATOMS);
for (j=FROZEN_ATOMS;jAtoms;j++)
for (k=0;k3;k++)
v[j][k]-=sum[k];
Sf=0.;
for (j=FROZEN_ATOMS;jAtoms;j++)
Sf+=v[j][0]*v[j][0]+v[j][1]*v[j][1]+v[j][2]*v[j][2];
/* Let's rescale the velocities */
Sf=sqrt(3*kb*temp*(Atoms-FROZEN_ATOMS)/Sf);
for (j=FROZEN_ATOMS;jAtoms;j++)
for (k=0;k3;k++)
v[j][k]*=Sf;
/*Now check it and be sure, that T appropriate and Vcm=0 */
Sf=0;
for (j=0;j3;j++)
sum[j]=0.;
for (j=FROZEN_ATOMS;jAtoms;j++){
Sf+=v[j][0]*v[j][0]+v[j][1]*v[j][1]+v[j][2]*v[j][2];
for (k=0;k3;k++){ sum[k]=sum[k]+v[j][k];}
}
fprintf(stdout,"\n Center_of_Mass
\n Vx=%lf Vy=%lf Vz=%lf ",sum[0],sum[1],sum[2]);
fprintf(stdout,"\n temp=%lf ", Sf/(3*(Atoms-FROZEN_ATOMS)*kb));
// for (j=0;jAtoms;j++)
//fprintf(stderr,"\n j=%ld %lf %lf %lf",j,v[j][0],v[j][1],v[j][2]);
}
////////////////////////////////////////////////
void bcc_001(void){
extern double r[Atoms][3];
extern double a0;
int i,j,k,pa;
double a02;//Initial coordinates
FILE *fd;
a02=a0/2;
pa=0;
for (k=0;kINIT_LAYERS;k++){
for (j=0;j size_x;j++){
for(i=0;isize_y;i++){
r[pa][0]=a0*i;
r[pa][1]=a0*j;
r[pa][2]=a02*k;
if(k%2!=0){
r[pa][0]+=a02;
r[pa][1]+=a02;}
pa++;
}
}
}
if((fd=fopen("r0.log","w"))!=NULL)
{
for (j=0;jAtoms;j++)
fprintf(fd,"\n %lf %lf %lf",
r[j][0],r[j][1],r[j][2]);
fclose(fd);
}
else
{fprintf(stderr,"\n I can't read file:\"
r0.log \"");exit(0);}
}
////////////////////////////
void bcc_011(void)
{
extern double r[Atoms][3];
extern double a0;
int i,j,k,n,index;
FILE *fd;
double E1[3][1],E2[3][1],E3[3][1],E4[3][1];
// The Unit vectors of the cell
E1[0][0]=a0*sqrt(2);E1[0][1]=E1[0][2]=0;
//E1=[sqrt(2);0;0];
E2[0][0]=E2[0][1]=a0/sqrt(2);E2[0][2]=0.;
//E2=[1./sqrt(2);1/sqrt(2);0];
E3[0][0]=E3[0][1]=0.;E3[0][2]=a0;
//E3=[0;0;1.0];
E4[0][0]=a0*sqrt(2)/2;E4[0][1]=0.;E4[0][2]=a0*0.5;
//E4=[sqrt(2)/2;0;1/2];
index=0;
for (k=0;ksize_y;k++){
for (j=0;jsize_z;j++){
for (i=0;isize_x;i++){
n=i+j*size_x+k*size_x*size_y;
r[index][2]=E1[1][0]*i+j*E2[1][0];
r[index][1]=E3[2][0]*k;
r[index][0]=E1[0][0]*i+E2[0][0];
// Set the center atoms for odd lines
if(j%2) r[index][0]+=E2[0][0];
index++;
}//i
}//j
}//k
/////////// ADD CENTRAL ATOM ////////////////
for (k=0;ksize_y;k++){
for (j=0;jsize_z;j++){
for (i=0;isize_x;i++){
r[index][2]=E4[1][0]*i+j*E2[1][0];
r[index][1]=E3[2][0]*k+E4[2][0];
r[index][0]=E1[0][0]*i+2*E4[0][0];
if(j%2) r[index][0]+=E4[0][0] ;
index++;
n++;
}//i
}//j
}//k
/////////////////////////////
if((fd=fopen("r0.log","w"))!=NULL)
{
for (j=0;jAtoms;j++)
fprintf(fd,"\n %lf %lf %lf",
r[j][0],r[j][1],r[j][2]);
fclose(fd);
}
else {fprintf(stderr,"\n I can't read file:\"
r0.log \"");exit(0);}
}
/////////////////////////////////////////////
void layer_colors(void)
{
extern double r[Atoms][3],Width,a0;
extern int Color[Atoms ]; double live[Atoms-FROZEN_ATOMS ][3], frozen[FROZEN_ATOMS ][3];
int i,j,cold=0,hot=0;
double h,barrier;
barrier=Width*(NUMBER_OF_FROZEN_LAYERS-1)+0.001;
for ( i=0;iAtoms;i++)
{
h=r[i][2]-barrier;
if(h=0)
{
frozen[cold][0]=r[i][0];frozen[cold][1]=r[i][1]; frozen[cold][2]=r[i][2];cold++;
//fprintf(stderr,"\n cold= %d Atoms=%d h=%g ", cold,Atoms,h);
}
else
{
live[hot][0]=r[i][0];live[hot][1]=r[i][1]; live[hot][2]=r[i][2];
hot++;
//fprintf(stderr,"\n hot= %d i=%d", hot,i);
}
}
// Chek post
if (cold-FROZEN_ATOMS)
{
fprintf(stderr, "That strange somehow !!!! cold != FROZEN_ATOMS");
fprintf(stderr, " \n \t cold =%d",cold);
fprintf(stderr, " \n \t FROZEN_ATOMS=%d",FROZEN_ATOMS);
exit(0);
}
if (hot+FROZEN_ATOMS-Atoms)
{fprintf(stderr, "That strange somehow !!!! hot != FROZEN_ATOMS-Atoms "); exit(0);}
//Rearrange the atoms:
for (i=0;icold;i++)
{ r[i][0]=frozen[i][0]; r[i][1]=frozen[i][1];r[i][2]=frozen[i][2];}
for (i=0;ihot;i++)
{ j=i+FROZEN_ATOMS;r[j][0]=live[i][0]; r[j][1]=live[i][1];r[j][2]=live[i][2];}
//Redefine Colors again:
for ( i=0;iAtoms;i++)
{
h=r[i][2]-Width/2.0;
if(h0) Color[i] =0;
else Color[i] =(int) h/Width+1;
}
//fprintf(stderr, "\n \n \t Bye !!!! ");
}
//////////////////////////////////////////
void check_initial_data()
{
switch(surface)
{
case 0: break;
case 1: break;
case 2: break;
default:
printf("This is an error, the parameter could be surface [0,3] only ");exit(0);
break;
}
/////////////////////////////////////////////
if(surface==1){
// Let's check the initial conditions:
if( 2*size_x*size_y*size_z-Atoms)
{
fprintf(stderr, "\n Attention : 2*size_x*size_y*size_z-Atoms
!=0 , \n Check the \" initdata.h
\" file ");
fprintf(stderr, "\n Write \"#define
Atoms %d \" ",2*size_x*size_y*size_z);
exit(0);}
if( 2*size_x*size_y*NUMBER_OF_FROZEN_LAYERS-FROZEN_ATOMS ) {
fprintf(stderr, "\n Attention : 2*size_x*size_y*NUMBER_OF_FROZEN_LAYERS-FROZEN_ATOMS
!= 0
, \n Check the \" initdata.h
\" file ");
fprintf(stderr, "\n Write \"
#define FROZEN_ATOMS %d \" ",
2*size_x*size_y*NUMBER_OF_FROZEN_LAYERS);
exit(0);
}
}
/////////////////////////
if(surface==2 surface==0 ){
if( size_x*size_y*size_z-Atoms)
{
fprintf(stderr, "\n Attention : size_x*size_y*size_z-Atoms
!=0 , \n Check the \" initdata.h
\" file ");
fprintf(stderr, "\n Write \"#define
Atoms %d \" ",size_x*size_y*size_z);
exit(0);}
if( size_x*size_y*NUMBER_OF_FROZEN_LAYERS-FROZEN_ATOMS )
{
fprintf(stderr, "\n Attention :
size_x*size_y*NUMBER_OF_FROZEN_LAYERS-FROZEN_ATOMS
!= 0 , \n Check the \" initdata.h
\" file ");
fprintf(stderr, "\n Write \"
#define FROZEN_ATOMS %d \" ",
size_x*size_y*NUMBER_OF_FROZEN_LAYERS);
exit(0);
}
}
}
///////////////////////////////
double lattice_parameter(double x)
{
int i;
double z;
// The Bulk simulations:
double lattice_set[10]=
{3.0482,3.0626,3.0706,3.073,3.089,3.092,3.095,3.096,3.096,3.18};
double temp_set[10]={400,1200, 1600,1700,2300,2400,2450,2480,2480,2500};
z= a00*exp(alpha_t*x);
if(x temp_set[0]){return(z);}
if(x temp_set[9]){return(z);}
if(x == temp_set[0]){return(lattice_set[0]);}
i=1;
while(xtemp_set[i])i++;
// The simplest linear interpolation:
z=(lattice_set[i]*(x-temp_set[i-1])+lattice_set[i-1]*(temp_set[i]-x))/
(temp_set[i]-temp_set[i-1]);
return(z);
}
////////////////////////// double scale_lattice() {
extern double a0;
extern double r[Atoms][3];
FILE *fd;
int i;
for (i=0;iAtoms;i++)
{r[i][0]*=a0;r[i][1]*=a0;r[i][2]*=a0;}
if((fd=fopen("r0.log","w"))!=NULL)
{
for( i=0;iAtoms;i++)
fprintf(fd," \n
%lf %lf %lf",r[i][0],r[i][1],r[i][2]);
fclose(fd);
}
}
***************''build_neigh_list.c.c'' *****************
Neighbor-list method - keep the list of the nearest neighbors
# include "constants.h"
# include "initdata.h"
# include "potential.h"
# include math.h
void build_neigh_list()
{
extern int nei[Atoms][Atoms];
extern double r[Atoms][3];
extern double lx,ly;
int i,j,k;
double dis[3]={0.,0.,0.};
double pr_dis;
for(i=0;iAtoms;i++)
for(j=0;jAtoms;j++)
nei[i][j]=0;
for(i=0;iAtoms;i++){
for(j=i+1;jAtoms;j++){
pr_dis=0.;
// x&y periodic boundary conditions
dis[0]=r[i][0]-r[j][0];
if(fabs(dis[0]) lx/2.0)
{
if(dis[0] 0.)dis[0]-=lx;
else dis[0]+=lx;
}
dis[1]=r[i][1]-r[j][1];
if(fabs(dis[1]) ly/2.0)
{
if(dis[1]0.)dis[1]-=ly;
else dis[1]+=ly;
}
dis[2]=r[i][2]-r[j][2];
pr_dis=sqrt(dis[0]*dis[0]+dis[1]*dis[1]+dis[2]*dis[2]);
//....................................................
if(pr_disneicutoff)
{
nei[i][0]++;
nei[i][nei[i][0]]=j;
}
} // for(i=0;iAtoms;i++)
} // for(j=i+1;jAtoms;j++)
}
***************''compute_f.c'' *****************
Calculation of force and potential energy per atom ( with the modified FS potential)
# include math.h
#include "initdata.h"
# include "func_cp.h"
//Modification of FS by Rebonato and et.al
#define Reb_FND 2.6319
#define Reb_K 3.3
#pragma_CRI inline density,density_dot,density_dot_dot,phi_pot,phi_pot_dot,
phi_pot_dot_dot,F_pot,F_pot_dot,F_pot_dot_dot
/* vectorization of compute_f,build_neigh_list , parallelization of compute_f */
//////////////////////////////////////////////////
void compute_f()
{
extern double r[Atoms][3],f[Atoms][3];
extern int regime,nei[Atoms][Atoms];
extern double pot_E,virsum,lx,ly;
extern double op;
int i,j,k,ll,mm;
double phi,dphi;
double dphi_ren,ddphi_ren;
double pr_dis,fr,fr1,fr2;
double rho,drho,ddrho;
double ui,dui,duj,ddui;
double dis[3],n[Atoms],Atoms];
/*Vectorization: A new part !!!! */
double nj[Atoms],f_tmp[Atoms][3],ni,fi0,fi1,fi2,pot_E_tmp,virsum_tmp;
pot_E=0.;
virsum=0.;
for(i=0;iAtoms;i++)
for(ll=0;ll3;ll++) f[i][ll]=0.0;
for(i=0;iAtoms;i++){
n[i]=0.0;
}
#pragma _CRI parallel shared(g11,g22,g23,n) defaults
for(j=0;jAtoms;j++) {
nj[j]=0.0;
}
/* The first loop: the background n(i) calculation */
#pragma _CRI taskloop
for( i=0;iAtoms;i++){
ni=0.0;
#pragma _CRI ivdep
for( k=1;k=nei[i][0];k++){
j=nei[i][k];
dis[0]=r[i][0]-r[j][0];
if (dis[0] lx/2.0 ) dis[0]-=lx;
if (dis[0] -lx/2.0) dis[0]+=lx;
dis[1]=r[i][1]-r[j][1];
if (dis[1] ly/2.0 ) dis[1]-=ly;
if (dis[1] -ly/2.0) dis[1]+=ly;
dis[2]=r[i][2]-r[j][2];
pr_dis=dis[0]*dis[0]+dis[1]*dis[1]+dis[2]*dis[2] ;
pr_dis=sqrt(pr_dis);
if(pr_discutoff)
{
rho=density(pr_dis);
//cj moved to the if regime block ! drho=density_dot(pr_dis);
//cj n[i]+=rho;
ni+=rho;
//cj n[j]+=rho;
nj[j]+=rho;
} // cj end of cutoff
} //cj end of k loop
// n[i]+=ni;
n[i]=ni;
} //cj end of i loop
#pragma _CRI guard
for(j=0;jAtoms;j++) {
n[j]+=nj[j];
if(regime==3){
g11[j]+=g11_tmp[j];
g22[j]+=g22_tmp[j];
g23[j]+=g23_tmp[j];
}
}
#pragma _CRI endguard
#pragma _CRI endparallel
//THE SECOND LOOP
#pragma _CRI parallel shared(g11,g22,g23,n) defaults
pot_E_tmp=0.0;
virsum_tmp=0.0;
for(j=0;jAtoms;j++)
for(ll=0;ll3;ll++) f_tmp[j][ll]=0.0;
#pragma _CRI taskloop
for( i=0;iAtoms;i++){
fi0=fi1=fi2=0.0;
ui=F_pot(n[i]);
dui=F_pot_dot(n[i]);
ddui=F_pot_dot_dot(n[i]);
pot_E_tmp+=ui;
#pragma _CRI ivdep
for( k=1;k=nei[i][0];k++){
j=nei[i][k];
dis[0]=r[i][0]-r[j][0];
if (dis[0] lx/2.0 ) dis[0]-=lx;
if (dis[0] -lx/2.0) dis[0]+=lx;
dis[1]=r[i][1]-r[j][1];
if (dis[1] ly/2.0 ) dis[1]-=ly;
if (dis[1] -ly/2.0) dis[1]+=ly;
dis[2]=r[i][2]-r[j][2];
pr_dis=dis[0]*dis[0]+dis[1]*dis[1]+dis[2]*dis[2] ;
pr_dis=sqrt(pr_dis);
if(pr_dis cutoff)
{
phi=phi_pot(pr_dis);
dphi=phi_pot_dot(pr_dis);
drho=density_dot(pr_dis);
if(regime==3)
{
dphi_ren=phi_pot_dot_ren(pr_dis);
ddphi_ren=phi_pot_dot_dot_ren(pr_dis);
ddrho=density_dot_dot(pr_dis);
}
pot_E_tmp+=phi;
duj=F_pot_dot(n[j]);
fr1=-dphi;
fr2=-drho*(dui+duj);
fr=fr1+fr2;
virsum_tmp+=fr*pr_dis;
for (ll=0;ll3;ll++)
{
//cj fi0,fi1,fi2 will be used instead of
f[i][ll]+=fr*dis[ll]/pr_dis;
f_tmp[j][ll]-=fr*dis[ll]/pr_dis;
}
fi0+=fr*dis[0]/pr_dis;
fi1+=fr*dis[1]/pr_dis;
fi2+=fr*dis[2]/pr_dis;
} //cutoff loop
} //nei loop
/*cj add fi0,fi1,fi2 to the force field */
f[i][0]=fi0;
f[i][1]=fi1;
f[i][2]=fi2;
} //i Atoms loop
#pragma _CRI guard
for(j=0;jAtoms;j++)
for(ll=0;ll3;ll++) f[j][ll]+=f_tmp[j][ll];
pot_E+=pot_E_tmp;
virsum+=virsum_tmp;
#pragma _CRI endguard
#pragma _CRI endparallel
}
/////////////////////////////////////////
double F_pot(double x)
{
return(A*sqrt(x));
}
///////////////////////////////////////////
double F_pot_dot(double x)
{
if(fabs(x)0.0000000000001)
x=0.0000000000001;
return(0.5*A/sqrt(x));
}
///////////////////////////////////////
double phi_pot(double x){
double y;
if(xr_cut) return(0);
y=(x-r_cut)*(x-r_cut)*(c1+c2*x+c3*x*x);
// FS+ Rebonato et al
if(xReb_FND)y+=Reb_K*(Reb_FND-x)*(Reb_FND-x)*(Reb_FND-x);
return(y);
}
/////////////////////////////////////////////// double phi_pot_dot(double x){
double y;
if(xr_cut)return(0);
y=2*(x-r_cut)*(c1+c2*x+c3*x*x)+(x-r_cut)*(x-r_cut)*(c2+2*c3*x);
// FS+ Rebonato et al
if(xReb_FND)y+=-3*Reb_K*(Reb_FND-x)*(Reb_FND-x);
return(y);
}
////////////////////////////////////////
double density(double pr_dis)
{
double rho;
if(pr_disFS_d)
return(0);
else
rho=(pr_dis-FS_d)*(pr_dis-FS_d);
return(rho);}
////////////////////////////////////////
double density_dot(double pr_dis)
{
double rho;
if(pr_disFS_d)
return(0);
else
rho=2*(pr_dis-FS_d);
return(rho);
}
*************** ``predictor_corrector.c '' *****************
Use the predictor corrector to solve the equations of motion
# include "initdata.h"
# include stdio.h
void corrector()
{
extern double r[Atoms][3],v[Atoms][3],f[Atoms][3];
extern double f1[Atoms][3],f2[Atoms][3],r1[Atoms][3],v1[Atoms][3];
extern double lx,ly;
extern int pbc[Atoms][2],regime;
double cr[]={3.0,10.,-1.};
double cv[]={7.,6.,-1.};
double div=24.;
int i;
for( i=FROZEN_ATOMS;iAtoms;i++){
r[i][0]=r1[i][0]+v1[i][0]*dt+(dt*dt/div)*(cr[0]*f[i][0]+cr[1]*f1[i][0]+
cr[2]*f2[i][0]);
r[i][1]=r1[i][1]+v1[i][1]*dt+(dt*dt/div)*(cr[0]*f[i][1]+cr[1]*f1[i][1]+
cr[2]*f2[i][1]);
r[i][2]=r1[i][2]+v1[i][2]*dt+(dt*dt/div)*(cr[0]*f[i][2]+cr[1]*f1[i][2]+
cr[2]*f2[i][2]);
v[i][0]=(r[i][0]-r1[i][0])/dt+(dt/div)*(cv[0]*f[i][0]+cv[1]*f1[i][0]+
cv[2]*f2[i][0]);
v[i][1]=(r[i][1]-r1[i][1])/dt+(dt/div)*(cv[0]*f[i][1]+cv[1]*f1[i][1]+
cv[2]*f2[i][1]);
v[i][2]=(r[i][2]-r1[i][2])/dt+(dt/div)*(cv[0]*f[i][2]+cv[1]*f1[i][2]+
cv[2]*f2[i][2]);
// Periodic Boundary Conditions for xy coordinates:
if(regime1)
{
if(r[i][0]lx){r[i][0]-=lx;pbc[i][0]++;}
if(r[i][1]ly){r[i][1]-=ly;pbc[i][1]++;}
if(r[i][0]0.){r[i][0]+=lx;pbc[i][0]-;}
if(r[i][1]0.){r[i][1]+=ly;pbc[i][1]-;}
}
else
{
if(r[i][0]lx){r[i][0]-=lx;}
if(r[i][1]ly){r[i][1]-=ly;}
if(r[i][0]0.){r[i][0]+=lx;}
if(r[i][1]0.){r[i][1]+=ly;}
}
}
}
//=================================//
void predictor()
{
extern double r[Atoms][3],v[Atoms][3],f[Atoms][3];
extern double f1[Atoms][3],f2[Atoms][3],r1[Atoms][3],v1[Atoms][3];
double cr[]={19.,-10.,3.};
double cv[]={27.,-22.,7.};
double div=24.;
int i,k;
for( i=FROZEN_ATOMS;iAtoms;i++){
for( k=0;k3;k++){
r1[i][k]=r[i][k];
v1[i][k]=v[i][k];
r[i][k]+=v[i][k]*dt+(dt*dt/div)*(cr[0]*f[i][k]+cr[1]*f1[i][k]+
cr[2]*f2[i][k]);
v[i][k]=(r[i][k]-r1[i][k])/dt+(dt/div)*(cv[0]*f[i][k]+cv[1]*f1[i][k]+
cv[2]*f2[i][k]);
f2[i][k]=f1[i][k];
f1[i][k]=f[i][k];
}
}
}
*************** ``evap_pto.c '' *****************
Evaluate the current pressure and temperature. Use the Noose-Hoover thermostat.
# include math.h
# include "initdata.h"
# include "constants.h"
void eval_PTO()
{
extern double ksi,lx,ly,lz,Q;
extern double f[Atoms][3],v[Atoms][3],r[Atoms][3];
extern double kin_E;
extern double temp,temp_current,press_current,virsum;
int i,k;
double vvsum,vol;
vol=lx*ly*lz;
vvsum=0.0;
for(i=FROZEN_ATOMS;iAtoms;i++)
vvsum+=(v[i][0]*v[i][0]+v[i][1]*v[i][1]+v[i][2]*v[i][2]);
kin_E=1.0*vvsum/2.;
temp_current=vvsum/(3.0*kb*(Atoms-FROZEN_ATOMS));
press_current=(temp_current*(Atoms-FROZEN_ATOMS)*kb+virsum/3.)/vol;
ksi+=dt*Q*((temp_current/temp)-1.0);
for(i=FROZEN_ATOMS;iAtoms;i++)
for(k=0;k3;k++)
f[i][k]=f[i][k]-Q*ksi*v[i][k];
}
*************** ``layers.c '' *****************
Evaluate the atom distribution (instant and average) between the layers. Calculate distance
between the neighboring layers.
# include "initdata.h"
# include "func_cp.h"
# include math.h
# include stdio.h
# include stdlib.h
void inter_layer_dist(int mode)
{
extern double r[Atoms][3],dis_bet_lay[NUMBER_OF_LAYERS-1],counter_dist;
extern int layers[NUMBER_OF_LAYERS][Atoms];
int i,k,p;
double z_prev, z_new;
FILE *fd;
if(mode==0)
{
//Initialization:
fprintf(stderr,"\n Inter Layer Distances initialization (:)");
counter_dist=0.;
for (i=0;iNUMBER_OF_LAYERS-1;i++)
dis_bet_lay[i]=0.;
}
//Start calculations
z_new=z_prev=0.;
for (i=0;iNUMBER_OF_LAYERS;i++)
{
z_new=0;
for (p=1;p=layers[i][0];p++)
{
k=layers[i][p];
z_new+=r[k][2];
}
if(layers[i][0]0)
{
z_new/=layers[i][0]; // Coordinate of c.m.c. of the layer
if(i0) dis_bet_lay[i-1]+=z_new-z_prev;
z_prev=z_new;
}
}
counter_dist++; // Number of the time measurements
if(mode==2)
{ //Result Output
if((fd=fopen("dis_bet_lay.txt","w"))!=NULL)
{
for (i=0;i(NUMBER_OF_LAYERS-1);i++)
{
if(counter_dist0)
fprintf(fd,"%d %lf \n",i+1,dis_bet_lay[i]/counter_dist);
}
fclose(fd);
}
}
}
//////////////////////////////////////
void inter_layer_dist_rw(char *fn,int mode)
{
extern double dis_bet_lay[NUMBER_OF_LAYERS-1],counter_dist;
int i;
FILE *fd;
if(mode==0)
{
////
if((fd=fopen(fn,"r"))!=NULL)
{
fscanf(fd,"\n%lf",&counter_dist);
for (i=0;iNUMBER_OF_LAYERS-1;i++)
fscanf(fd,"\n%lf", &dis_bet_lay[i]);
fclose(fd);
}
else inter_layer_dist(0);
////
}
if(mode==1)
{
////////////
if((fd=fopen(fn,"w"))!=NULL)
{
fprintf(fd,"\n %lf",counter_dist);
for (i=0;iNUMBER_OF_LAYERS-1;i++)
fprintf(fd,"\n%lf", dis_bet_lay[i]);
fclose(fd);
}
///////////
}
}
//=============================================//
void atom_distribution_in_layers (int mode)
{
extern double r[Atoms][3],Width;
extern int layers[NUMBER_OF_LAYERS][Atoms];
//[][]-Number of a layer, number of Atoms [0][],[][i0] indexes
extern double a0;
extern int atoms_distribution[NUMBER_OF_LAYERS];
extern long counter_of_atoms_distribution,step;
double Z0;
double dz;
int i,j,index;
FILE * fd;
char fn[80];
if(mode==0)
{
fprintf(stderr,"\n Atoms distribution in the layers initialization (:) ");
counter_of_atoms_distribution=0;
for (i=0;iNUMBER_OF_LAYERS;i++)
atoms_distribution[i] =0;
}
for (i=0;iNUMBER_OF_LAYERS;i++)
for (j=0;jAtoms;j++)
layers[i][j]=0;
// Collect the Statistics
Z0=NUMBER_OF_FROZEN_LAYERS *Width-Width/2;
for (i=FROZEN_ATOMS;iAtoms;i++)
{
dz=r[i][2]-Z0;
if(dz0)
{
if(r[i][2](NUMBER_OF_FROZEN_LAYERS -1)*Width )
{
fprintf(stderr,"\n layers.c: Very Negative Base Layer
i=%d Z0=%lf r(i,3)=%lf diff=%lf ",i,Z0,r[i][2],r[i][2]-Z0);
exit(0);
}
dz=-dz;
}
//////////
if(fmod(dz,Width)==0)
{
index=(int)(dz/Width)-1;
if(indexNUMBER_OF_LAYERS){layers[index][0]+=1; layers[index][layers[index][0]]=i;}
else {fprintf(stderr,"\n layers.c :Super adlayer %d",index);exit(0);}
}
else
{index=(int)(dz/Width);
if(indexNUMBER_OF_LAYERS)
{
layers[index][0]+=1;
layers[index][layers[index][0]]=i;
}
else {
//Adlayer
index=NUMBER_OF_LAYERS-1;
layers[index][0]+=1;
layers[index][layers[index][0]]=i;
}
}
}
for (i=0;iNUMBER_OF_LAYERS;i++)
{
atoms_distribution[i] +=layers[i][0];
sprintf(fn,"layer_%d.txt",i);
if((fd=fopen(fn,"a+"))!=NULL){
fprintf(fd,"\n%ld %d", step,layers[i][0]);
fclose(fd); }
}
counter_of_atoms_distribution++;
}
//=========================================================//
void atom_distribution_in_layers_rw (char *fn, int mode)
{
extern int atoms_distribution[NUMBER_OF_LAYERS];
extern long counter_of_atoms_distribution;
int i;
FILE * fd;
if(mode==0)
{
if((fd=fopen(fn,"w"))!=NULL)
{
///////////////
fprintf(fd,"\n%ld",counter_of_atoms_distribution);
for (i=0;iNUMBER_OF_LAYERS;i++)
fprintf(fd,"\n%d",atoms_distribution[i] );
fclose(fd);
//////////////////
}
}
if(mode==1)
{
if((fd=fopen(fn,"r"))!=NULL)
{
fscanf(fd,"\n%ld",&counter_of_atoms_distribution);
for (i=0;iNUMBER_OF_LAYERS;i++)
fscanf(fd,"\n%d", &atoms_distribution[i] );
fclose(fd);
}
else atom_distribution_in_layers (0);
}
}
*************** ``density_profile.c'' *****************
Evaluate the atomic density along the z-direction. Use the Chen et al. method to facilitate
the obtained data.
# include stdio.h
# include math.h
# include stdlib.h
# include "func_cp.h"
# include "potential.h"
# include "initdata.h"
# include "constants.h"
void chen_density_profile(int mode)
{
extern double r[Atoms][3],Width;
extern double dens_ch[Z_RES],chen_counter;
double z_in,z_fin,delta,z,f,sigma_ch;
int i,j;
FILE *fd;
z_in=NUMBER_OF_FROZEN_LAYERS*Width;
z_fin=z_in+(NUMBER_OF_LAYERS+4)*Width;
delta=(z_fin-z_in)/Z_RES;
sigma_ch=0.1*Width;
if(mode==0)
{
//Initialization:
fprintf(stderr,"\n Chen density profile initialization (:) ");
chen_counter=0;
for (i=0;iZ_RES;i++)
dens_ch[i]=0.;
}
// Let's start calculations
for (i=0;iZ_RES;i++)
{
z=z_in+i*delta;
for (j=FROZEN_ATOMS;jAtoms;j++)
{
f=(z-r[j][2])*(z-r[j][2])/(2*sigma_ch*sigma_ch);
dens_ch[i]+=exp(-f)/(sqrt(2*pi)*sigma_ch);
}
}
chen_counter++; // Number of the measurements
// Let's print the results
if(mode==2)
{
if((fd=fopen("chen_density.txt","w"))!=NULL)
{
for (i=0;iZ_RES;i++)
fprintf(fd,"%lf %lf \n",z_in+i*delta,dens_ch[i]/chen_counter);
fclose(fd);
}
}
}
/////////////////////////////
void chen_density_profile_rw(char *fn, int mode)
{
extern double dens_ch[Z_RES],chen_counter;
int i;
FILE *fd;
if(mode==1)
{
// BackUp
if((fd=fopen(fn,"w"))!=NULL)
{
fprintf(fd,"%lf",chen_counter);
for (i=0;iZ_RES;i++) fprintf(fd,"%lf ",dens_ch[i]);
fclose(fd);
}
}
if(mode==0)
{
//Read
if((fd=fopen(fn,"r"))!=NULL)
{
fscanf(fd,"%lf",&chen_counter);
for (i=0;iZ_RES;i++) fscanf(fd,"%lf ",&dens_ch[i]);
fclose(fd);
}
else chen_density_profile(0);
}
}
*************** ``temp_prof.c '' *****************
Find the temperature profile along the z-directions.
# include "initdata.h"
# include "constants.h"
# include "func_cp.h"
# include "potential.h"
# include math.h
# include stdio.h
void temperature_profile (int mode)
{
extern double v[Atoms][3],temp_profile[NUMBER_OF_LAYERS],counter_temp;
extern int layers[NUMBER_OF_LAYERS][Atoms];
int i,j,k;
double s;
FILE *fd;
if(mode==0)
{ //Initialization:
fprintf(stderr,"\n Temperature profile initialization (:) ");
counter_temp=0.;
for (i=0;iNUMBER_OF_LAYERS;i++)
temp_profile[i]=0.;
}
//Start calculations
for (i=0;iNUMBER_OF_LAYERS;i++)
{
s=0.;
for (j=1;j=layers[i][0];j++)
{
k=layers[i][j];
s+=(v[k][0]*v[k][0]+v[k][1]*v[k][1]+v[k][2]*v[k][2]);
}
if(layers[i][0]10)
{
s/=kb*3*layers[i][0];
temp_profile[i]+=s;
}
}
counter_temp++; // Number of the time measurements
if(mode==2)
{ //Result Output
if((fd=fopen("temp_profile.txt","w"))!=NULL)
{
for (i=0;iNUMBER_OF_LAYERS;i++)
if(counter_temp0) fprintf(fd,"%d %lf \n",i+1,temp_profile[i]/counter_temp);
// for (i=0;iNUMBER_OF_LAYERS;i++)
// fprintf(stderr," \n %lf %lf ",counter_temp,temp_profile[i]/counter_temp);
fclose(fd);
}
}
}
//======================//
void temperature_profile_rw(char *fn,int mode)
{
extern double temp_profile[NUMBER_OF_LAYERS],counter_temp;
int i;
FILE *fd;
if(mode==0)
{
if((fd=fopen(fn,"r"))!=NULL)
{
for (i=0;iNUMBER_OF_LAYERS;i++)
fscanf(fd," \n %lf ",&temp_profile[i]);
fscanf(fd," \n %lf ",&counter_temp);
fclose(fd);
}
else temperature_profile( 0 ) ; // initialization
}
if(mode==1)
{
if((fd=fopen(fn,"w"))!=NULL)
{ for (i=0;iNUMBER_OF_LAYERS;i++)
fprintf(fd," \n %lf ",temp_profile[i]);
fprintf(fd," \n %lf ",counter_temp);
fclose(fd);
}
}
}
*************** ``order_parameter.c '' *****************
Calculate the structure order parameters along the x,y, and z-directions.
# include "initdata.h"
# include "constants.h"
# include "func_cp.h"
# include "potential.h"
# include math.h
# include stdio.h
void order_parameters(int mode)
{
extern double r[Atoms][3];
extern int layers[NUMBER_OF_LAYERS][Atoms];
// The order parameter is defined in each layer and in various directions
extern double S[NUMBER_OF_LAYERS][NUMBER_OF_DIRECTIONS][2] ;
extern double a0;
double k[NUMBER_OF_DIRECTIONS][3];
double sx,sy,sz,dx,dy,S_c,S_s;
int index,p,i,j,Ns,parity;
double or1,or2,or3;
FILE *fd;
char fn[50];
// Directions [001]
if(surface==0)
{
k[0][0]=4*pi/a0; k[0][1]=0.; k[0][2]=0.;
k[1][0]=0.; k[1][1]=4*pi/a0; k[1][2]=0.;
k[2][0]=0.; k[2][1]=0.; k[2][2]=4*pi/a0;
}
if(surface==1)
{
k[0][0]=2*sqrt(2)*pi/a0; k[0][1]=0.; k[0][2]=0;
k[1][0]=0; k[1][1]=4*pi/a0; k[1][2]=0.;
k[2][0]=0; k[2][1]=0; k[2][2]=2*sqrt(2)*pi/a0;
}
if(surface==2)
{
k[0][0]=4*pi/(a0*sqrt(2.)); k[0][1]=0.; k[0][2]=0.;
k[1][0]=0; k[1][1]=4*pi/(2.*a0*sqrt(3./2.)); k[1][2]=0.;
k[2][0]=0; k[2][1]=0; k[2][2]=2*pi/(a0/(2*sqrt(3.)));
}
if(mode ==0)
{
// Initiallization of the arrays
fprintf(stderr,"\n Order parameter profile initialization (:) ");
for (i=0;iNUMBER_OF_LAYERS;i++)
for (j=0;jNUMBER_OF_DIRECTIONS;j++)
{ S[i][j][0]=0.; S[i][j][1]=0.; }
}
for (p=0;p3;p++)//3D space
{//p
for (i=0;iNUMBER_OF_LAYERS;i++)
{//Initialization of the arrays
Ns=0;
S_c=S_s=0.0;
for (j=1;j=layers[i][0];j++)
{
index=layers[i][j];
//Take all the particles
sz=k[p][0]*r[index][0]+k[p][1]*r[index][1]+k[p][2]*r[index][2];
S_c+=cos(sz);
S_s+=sin(sz);
Ns++;
} // end of the j cycle
if(Ns)sx=(S_c * S_c+S_s * S_s)/ (Ns*Ns);
if(layers[i][0]20)
{
S[i][p][0]+=sx; S[i][p][1]++;
}
} // end of the i Layer
} // 3D space
if(mode==2)
{
if((fd=fopen("order_prof.txt","w"))!=NULL){
for (i=0;iNUMBER_OF_LAYERS;i++)
{
or1=or2=or3=0.;
if(S[i][0][1]0)or1=S[i][0][0]/S[i][0][1];
if(S[i][1][1]0)or2=S[i][1][0]/S[i][1][1];
if(S[i][2][1]0)or3=S[i][2][0]/S[i][2][1];
fprintf(fd,"\n %d %lf %lf %lf",i+1,or1,or2,or3);
// Start calculate from 1
}
fclose(fd);
}
else{fprintf(stderr,"\n Can't open the %s file ",fn);exit(0);}
}//if(mode==2)
}
//====================================================//
void order_parameters_rw(char *fn, int mode )
{
extern double S[NUMBER_OF_LAYERS][NUMBER_OF_DIRECTIONS][2] ;
int i;
FILE *fd;
if(mode==0)
{
//////////////////////////////////////
if((fd=fopen(fn,"r"))!=NULL)
{
for (i=0;iNUMBER_OF_LAYERS;i++)
fscanf(fd," \n %lf %lf %lf %lf %lf %lf",
&S[i][0][0],&S[i][1][0],&S[i][2][0],&S[i][0][1],&S[i][1][1],&S[i][2][1]);
fclose(fd);
}
else order_parameters(0);
//////////////////////////////////////
}
if(mode==1)
{
//////////////////////////////////////
if((fd=fopen(fn,"w"))!=NULL)
{
for (i=0;iNUMBER_OF_LAYERS;i++)
fprintf(fd,"\n %lf %lf %lf %lf %lf %lf",
S[i][0][0],S[i][1][0],S[i][2][0],S[i][0][1],S[i][1][1],S[i][2][1]);
fclose(fd);
}
//////////////////////////////////////
}
}
*************** `` pair_correlation.c '' *****************
Evaluate the atomic 2D radial distribution function
# include "initdata.h"
# include "constants.h"
# include "func_cp.h"
# include "potential.h"
# include math.h
# include stdio.h
# include stdlib.h
void pair_correlation(int mode)
{
extern double r[Atoms][3],Profile[NUMBER_OF_LAYERS][PAIR_STEP];
extern double lx,ly,width_rdf, normFac,half_box;
extern int layers[NUMBER_OF_LAYERS][Atoms];
extern long pair_corr_count;
extern double a0;
int i, j, k, index;
int kk,jj;
double dist,s;
double dis[2];
char fn[30];
FILE * fd;
if(mode ==0)
{
if (lxly) half_box=ly/2;
else half_box=lx/2;
width_rdf=half_box/PAIR_STEP;
normFac=lx*ly/(2*pi*width_rdf*width_rdf);
// Initiallization of the arrays
fprintf(stderr,"\n Pair Correlation initialization (:) ");
for (i=0;iNUMBER_OF_LAYERS;i++)
for (j=0;jPAIR_STEP;j++)
Profile[i][j]=0.;
pair_corr_count=0;
}
// Designing the G(r) function
for (i=0;iNUMBER_OF_LAYERS;i++)
{
for (kk=1;kklayers[i][0];kk++) { // layers[i][0] is instant number of atoms in i layer
for (jj=kk+1;jj=layers[i][0];jj++) {
//Periodic Boundary Conditions for the xy (parallel to the surafce) components
dist=0.;
k=layers[i][kk];
j=layers[i][jj];
dis[0]=r[j][0]-r[k][0];
if(fabs(dis[0]) lx/2.0)
{
if(dis[0]0.)dis[0]-=lx;
else dis[0]+=lx;
}
dis[1]=r[j][1]-r[k][1];
if(fabs(dis[1]) ly/2.0)
{
if(dis[1]0.)dis[1]-=ly;
else dis[1]+=ly;
}
// End of the Periodic Boundary Conditions for the xy (paralle to the surafce) components
dist=dis[0]*dis[0]+dis[1]*dis[1];
dist=sqrt(dist);
// Check what you still inside the box (just to be on the safe side)
index=(int)(dist/width_rdf);
if ( index PAIR_STEP ) {
if(layers[i][0]0){
Profile[i][index]+=2./(layers[i][0]*layers[i][0]); // Rij & Rji contributions
}
} // j j
} // kk
} // i: NUMBER_OF_LAYERS
pair_corr_count++;
} //==================
if(mode==2){
/*Final Output for each layer*/
for (i=0;iNUMBER_OF_LAYERS;i++)
{
sprintf(fn,"pair_func_%d.txt",i);
if((fd=fopen(fn,"w"))!=NULL)
{
for(j=0;jPAIR_STEP;j++)
{
if(pair_corr_count0){
fprintf(fd,"%lf %lf \n",(j+0.5)*width_rdf, normFac*Profile[i][j]/((j+0.5)*pair_corr_count));
}
else
{fprintf(fd,"%lf %lf \n",j*width_rdf,0.); }
}
fclose(fd); }
} // end for
}
}
//===============================//
void pair_correlation_rw(char *fn, int mode )
{
extern double Profile[NUMBER_OF_LAYERS][PAIR_STEP];
extern long pair_corr_count;
int i,j;
FILE *fd;
if(mode==0)
{
if((fd=fopen(fn,"r"))!=NULL)
{
for (i=0;iNUMBER_OF_LAYERS;i++)
for (j=0;jPAIR_STEP;j++)
fscanf(fd," \n %lf ",&Profile[i][j]);
fscanf(fd," \n %ld ",&pair_corr_count);
fclose(fd);
}
else pair_correlation(0);
}
if(mode==1)
{
if((fd=fopen(fn,"w"))!=NULL)
{
for (i=0;iNUMBER_OF_LAYERS;i++)
for (j=0;jPAIR_STEP;j++)
fprintf(fd," \n %lf ",Profile[i][j]);
fprintf(fd," \n %ld ",pair_corr_count);
fclose(fd);
}
}
}
*************** ``diffusion.c '' *****************
Calculate in-plane the out-plane diffusion coefficients
diff_counter - number of diffusion(1) calls
NUMBER_OF_CALLS - number of calls when a new t0 is taken
interval - number of steps betwee two succesive calls
T0max - maximal number of time origins in the ``initdata.h''
TMAX - maximal number of sampling in the ``initdata.h''
t0max - maximal number of time origins
tmax - maximal number of sampling
t0 - number of time origins
# include "initdata.h"
# include "constants.h"
# include "func_cp.h"
# include "potential.h"
# include math.h
# include stdio.h
void diffusion(int mode)
{
extern double r[Atoms][3],lx,ly;
extern int layers[NUMBER_OF_LAYERS][Atoms];
extern long num_of_meas,diff_counter,t0max,tmax;
extern int ntime[TMAX],pbc[Atoms][2],time0[T0max],t0,interval;
extern double r0[Atoms][3][T0max];
extern double r2f[TMAX][3][NUMBER_OF_LAYERS];
long delta,tt0;
int i,j,k,n,p;
double dtime, norm;
FILE *fr;
char fn[30];
if(mode ==0)
{
fprintf(stderr,"\n Diffusion initialization (:) ");
tmax=t0+num_of_meas;
if (tmax TMAX)
{
if((fr=fopen("README","w"))!=NULL)
{
fprintf(fr,"\n tmax=%ld TMAX=%ld",tmax,TMAX);
fprintf(fr,"\n Only the shortest time TMAX=%ld is taken ",TMAX);
fclose(fr);
tmax=TMAX;
}
}
t0max=(long)1.*tmax/NUMBER_OF_CALLS;
if (t0max T0max) t0max = T0max;
for(i=t0;it0max;i++)
time0[i]=0;//The current time origin
for(i=t0;itmax;i++) {
ntime[i]=0; // delta table
for(j=0;jNUMBER_OF_LAYERS;j++) {
for(k=0;k3;k++)
r2f[i][k][j]=0.; }
}
for (i=FROZEN_ATOMS;iAtoms;i++)
for (k=t0;kt0max;k++)
for (j=0;j3;j++)
r0[i][j][k]=0.;//Starting Points
} // End of Start Initialization
// Start of Sampling
if(diff_counter%NUMBER_OF_CALLS==0)
{ //Let's take a new time origin
tt0=t0-(int)(t0/t0max)*t0max;//Check if number of time origins less than tmax, if not remove the first meas.
++t0; //Number of time origins are taken
time0[tt0]=diff_counter; // Store the time at t=0;
for (i=FROZEN_ATOMS;iAtoms;i++)
{
//For Future Parallelization
r0[i][0][tt0]=r[i][0]+pbc[i][0]*lx; //Starting Points
r0[i][1][tt0]=r[i][1]+pbc[i][1]*ly; //Starting Points
r0[i][2][tt0]=r[i][2]; //Starting Points
} // for(i=FROZEN_ATOMS;iAtoms;i++)
} // if(diff_counter%NUMBER_OF_CALLS==0)
n=(t0t0max)?t0:t0max;
for(k=0;kn;k++) //Update the r2f now
{
delta=diff_counter-time0[k]; //Actual time minus time origin
if(deltatmax) // Watch out the matrix boundaries
{
ntime[delta]++;
// norm=1.0/(Atoms-FROZEN_ATOMS);// Divide on the number of the particles in the layer
for (j=0;jNUMBER_OF_LAYERS;j++)
{
for (p=1;p=layers[j][0];p++)
{
i=layers[j][p];
if(layers[j][0]) norm=1.0/layers[j][0];// Divide on the number of the particles in the layer
else norm=0.0;
r2f[delta][0][j]+=(r[i][0]+pbc[i][0]*lx-r0[i][0][k])*(r[i][0]+pbc[i][0]*lx
- r0[i][0][k])*norm;
r2f[delta][1][j]+=(r[i][1]+pbc[i][1]*ly-r0[i][1][k])*(r[i][1]+pbc[i][1]*ly
- r0 [i][1][k])*norm;
r2f[delta][2][j]+=(r[i][2]-r0[i][2][k])*(r[i][2]-r0[i][2][k])*norm;
} // for (p=1;p=layers[j][0];p++)
} // for (j=0;jNUMBER_OF_LAYERS;j++)
}// if(deltaTMAX)
} // for(k=0;kn;k++)
diff_counter++; // Number of calls
// End of the Sampling
if(mode==2)
{ // Output of the Results
dtime=2*0.721*dt*interval;
for (k=0;kNUMBER_OF_LAYERS;k++){
sprintf(fn,"r2f%d.txt",k);
if((fr=fopen(fn,"w"))!=NULL){///----
for(i=0;itmax;i++)
{
if(ntime[i]0)
{
norm=1./ntime[i];
fprintf(fr,"%lf %lf %lf %lf \n",
dtime*(i+0.5),r2f[i][0][k]*norm,r2f[i][1][k]*norm,r2f[i][2][k]*norm);
} //(ntime[i]0)
} // for(i=0;iTMAX;i++)
fclose(fr); }
else{fprintf(stderr,"\n\n\n Can't open: r2f%d.txt",k);}
}
}//End of Output of the Results }
//====================//
void diffusion_rw(char *fn,int mode)
{
extern long diff_counter,t0max,tmax;
extern int ntime[TMAX],time0[T0max],t0;
extern double r0[Atoms][3][T0max];
extern double r2f[TMAX][3][NUMBER_OF_LAYERS];
FILE *fd;
int i,j;
if (mode==0)
{
if((fd=fopen(fn,"r"))!=NULL)
{
fscanf(fd," %ld %ld %ld ",&diff_counter,&t0max,&tmax);
fscanf(fd," \n %d ",&t0);
for(i=0;itmax;i++)
fscanf(fd," \n %d ",&ntime[i]);
for(i=0;it0max;i++)
fscanf(fd," \n %d ",&time0[i]);
for(i=0;itmax;i++)
for(j=0;jNUMBER_OF_LAYERS;j++)
fscanf(fd," \n %lf %lf %lf ",&r2f[i][0][j],&r2f[i][1][j],&r2f[i][2][j]);
for (i=FROZEN_ATOMS;iAtoms;i++)
for (j=0;jt0max;j++)
fscanf(fd," \n %lf %lf %lf ",&r0[i][0][j],&r0[i][1][j],&r0[i][2][j]);
fclose(fd);
}
else {diffusion(0);system("rm -f pbc.back;");}
}
if (mode==1)
{
if((fd=fopen(fn,"w"))!=NULL)
{
fprintf(fd," %ld %ld %ld ",diff_counter,t0max,tmax);
fprintf(fd," \n %d ",t0);
for(i=0;itmax;i++)
fprintf(fd," \n %d ",ntime[i]);
for(i=0;it0max;i++)
fprintf(fd," \n %d ",time0[i]);
for(i=0;itmax;i++)
for(j=0;jNUMBER_OF_LAYERS;j++)
fprintf(fd," \n %lf %lf %lf ",r2f[i][0][j],r2f[i][1][j],r2f[i][2][j]);
for (i=FROZEN_ATOMS;iAtoms;i++)
for (j=0;jt0max;j++)
fprintf(fd," \n %lf %lf %lf ",r0[i][0][j],r0[i][1][j],r0[i][2][j]);
fclose(fd);
} }
}
*************** ``read_save.c '' *****************
The various input-output functions:
Aviz, mixing, read_vel_from_file(char *fn), read_coord_from_file(char *fn )
save_to_file(long no,int d), save_inf(long step,long counter,float stam_real), get_inf(void)
pbc_recover(void), pbc_init(void)
#include stdio.h
#include stdlib.h
#include string.h
# include "initdata.h"
# include "constants.h"
# include "func_cp.h"
void aviz(char *fn)
{
extern double r[Atoms][3];
extern int Color[Atoms ];
FILE * fd;
int j,i;
char *color[]=
{"a1","a2","a3","a4","a5",
"a6","a7","a8","a9","b1",
"b2","b3","b4","b5","b6",
"b7","b8","b9","c1","c2",
"c3","c4","c5","c6","c7",
"c8","c9","d1","d2","d3"} ;
if(surface==0) {
if((fd=fopen(fn,"w"))!=NULL)
{
fprintf(fd," %d ", Atoms);
fprintf(fd,"\n ### Surface simulations ");
for (j=0;jFROZEN_ATOMS;j++)
fprintf(fd,"\n Cu %lf %lf %lf ", r[j][0],r[j][1],r[j][2]);
for(i=0;iINIT_LAYERS-NUMBER_OF_FROZEN_LAYERS;i++)
for (j=FROZEN_ATOMS+i*size_x*size_y;j FROZEN_ATOMS+(i+1)*size_x*size_y; j++)
fprintf(fd,"\n %s %lf %lf %lf ",color[i],r[j][0],r[j][1],r[j][2]);
fclose(fd);
}
else { fprintf(stderr,"I can't open the file %s !!!",fn); exit(0);}
} if(surface==1)
{
if((fd=fopen(fn,"w"))!=NULL)
{
fprintf(fd," %d ", Atoms);
fprintf(fd,"\n ### Surface simulations [011] ");
for (j=0;jAtoms;j++)
if ( Color[j] NUMBER_OF_FROZEN_LAYERS)fprintf(fd,"\n Cu %lf %lf %lf ", r[j][0],r[j][1],r[j][2]);
else fprintf(fd,"\n %s %lf %lf %lf ",color[Color[j]],r[j][0],r[j][1],r[j][2]);
fclose(fd);
}
else { fprintf(stderr,"I can't open the file %s !!!",fn); exit(0);}
}
if(surface==2) {
if((fd=fopen(fn,"w"))!=NULL)
{
fprintf(fd," %d ", Atoms);
fprintf(fd,"\n ### Surface simulations [111] ");
for (j=0;jAtoms;j++)
if ( Color[j] NUMBER_OF_FROZEN_LAYERS)fprintf(fd,"\n Fe %lf %lf %lf ", r[j][0],r[j][1],r[j][2]);
else fprintf(fd,"\n %s %lf %lf %lf ",color[Color[j]],r[j][0],r[j][1],r[j][2]);
fclose(fd);
}
else { fprintf(stderr,"I can't open the file %s !!!",fn); exit(0);}
}
}
//===============================//
void read_vel_from_file(char *fn)
{
extern double v[Atoms][3];
int i;
FILE *fp;
fprintf(stderr,"\n %s",fn);
if((fp=fopen(fn,"r"))!=NULL)
{
for(i=0;iAtoms;i++)
fscanf(fp,"%lf %lf %lf \n",&v[i][0],&v[i][1],&v[i][2]);
//printf("N=%d %g %g %g \n",i,v[i][0],v[i][1],v[i][2]);
fclose(fp);
}
else
{fprintf(stderr,"\n I can't read file: %s",fn);exit(0);}
}
//======================================//
void read_coord_from_file(char *fn)
{
extern double r[Atoms][3];
int i;
FILE *fd;
fprintf(stderr,"\n %s",fn);
if((fd=fopen(fn,"r"))!=NULL){
for(i=0;iAtoms;i++)
fscanf(fd,"%lf %lf %lf \n",&r[i][0],&r[i][1],&r[i][2]);
//fprintf(stderr,"N=%d %g %g %g \n",i,r[i][0],r[i][1],r[i][2]);
fclose(fd); }
else
{fprintf(stderr,"\n I can't read file: %s",fn);exit(0);}
}
//==========================//
void energy_temp_pressure()
{
extern long step;
extern double pot_E,kin_E,temp;
extern double temp_current,press_current;
FILE *fd;
if((fd=fopen("infoEPT.txt","a+"))!=NULL)
{
fprintf(fd,"\n %ld %lf %lf %lf %lf",step,pot_E,kin_E,temp_current,press_current);
fclose(fd);
}
}
//==========================================//
void save_configuration(int mode)
{
extern double r[Atoms][3],v[Atoms][3];
FILE *fd;
int i;
if(mode)
{
if((fd=fopen("r.log","w"))!=NULL)
{
for( i=0;iAtoms;i++)
fprintf(fd," \n %lf %lf %lf",r[i][0],r[i][1],r[i][2]);
fclose(fd);
}
if((fd=fopen("v.log","w"))!=NULL)
{
for( i=0;iAtoms;i++)
fprintf(fd," \n %lf %lf %lf",v[i][0],v[i][1],v[i][2]);
fclose(fd);
}
}
}
//===================================//
void save_inform(char *fn)
{
extern int regime,interval,Limit_aviz_pict;
extern long step,step0,num_eq_steps,num_of_meas;
extern double temp,Q;
FILE *outd;
float percent;
// Calculate percentage of the accomplished job
if(regime2)
{percent=100.0*(step-step0)/num_eq_steps;}
else
{percent=100.0*(step-step0)/(num_of_meas*interval);}
if((outd=fopen(fn,"w"))!=NULL)
{
fprintf(outd,"\n %d",regime);
fprintf(outd,"\n %4.0f",temp);
fprintf(outd,"\n %ld",num_eq_steps);
fprintf(outd,"\n %d",interval);
fprintf(outd,"\n %ld",num_of_meas);
fprintf(outd,"\n %ld",step);
fprintf(outd,"\n %d",Limit_aviz_pict);
fprintf(outd,"\n %3.1f",percent);
fprintf(outd,"\n %f",Q);
fprintf(outd,"\n %d",Atoms);
fprintf(outd,"\n %d",(int)size_x*size_y);
fprintf(outd,"\n %d",NUMBER_OF_FROZEN_LAYERS);
fprintf(outd,"\n %d",INIT_LAYERS-NUMBER_OF_FROZEN_LAYERS);
switch (surface)
{
case 0:
fprintf(outd,"\n [001]");
break;
case 1:
fprintf(outd,"\n [011]");
break;
case 2:
fprintf(outd,"\n [111]");
break;
default:
printf("Too Strange, too much different surface=%d \n", surface);
exit(0);
break;
}
fclose(outd);
comments(fn);
}
else
{
fprintf(stderr,"\n Can't open the file:\"%s\"",fn);
exit(0);
}
}
//===========================//
long get_inform(char *fn)
{
// Get all necessary information
extern int regime,interval,Limit_aviz_pict,step_Aviz;
extern long step,step0,num_eq_steps,num_of_meas;
extern double temp,Q;
float percent;
FILE *outd;
if((outd=fopen(fn,"r"))!=NULL)
{
fscanf(outd,"\n %d",®ime);
fscanf(outd,"\n %lf",&temp);
fscanf(outd,"\n %ld",&num_eq_steps);
fscanf(outd,"\n %d",&interval);
fscanf(outd,"\n %ld",&num_of_meas);
fscanf(outd,"\n %ld",&step0);
fscanf(outd,"\n %d",&Limit_aviz_pict);
fscanf(outd,"\n %g",&percent);
fscanf(outd,"\n %lf",&Q);
fclose(outd);
if (regime) fprintf(stderr,"\n Regime: %d - measurements",regime);
else fprintf(stderr,"\n Regime: %d - equilibration",regime);
fprintf(stderr,"\n\n Temperature %4.0f K",temp);
fprintf(stderr,"\n Number of steps in approaching to the equilibrium %ld",num_eq_steps);
fprintf(stderr,"\n Number of steps in the measurement regime %d \n ",interval*num_of_meas);
fprintf(stderr,"\n\t Number of steps between two measurements %d",interval);
fprintf(stderr,"\n\t Number of measurements %ld",num_of_meas);
fprintf(stderr,"\n \t The last step was = %ld \n ",step0);
fprintf(stderr,"\n \t The task have been complited on = %3.0f percent",percent);
fprintf(stderr,"\n \n \t Layers + addlayer = %d + %d",NUMBER_OF_LAYERS-1,1);
fprintf(stderr,"\n \t Frozen Atoms =%d",FROZEN_ATOMS);
fprintf(stderr,"\n \t Frozen Layers = %d",NUMBER_OF_FROZEN_LAYERS );
fprintf(stderr,"\n \n \t Total number of steps =%d",num_eq_steps+num_of_meas*interval);
fprintf(stderr,"\n \t Real time = %lf psec",(num_eq_steps+num_of_meas*interval)*dt*0.721);
step_Aviz=(int)(1.0*num_eq_steps+num_of_meas*interval)/Limit_aviz_pict;
}
else
{
fprintf(stderr,"\n Can't open the file:\"%s\" , I try to guess ",fn);
if((outd=fopen(fn,"w"))!=NULL)
{
fprintf(outd,"\n %d",0);
fprintf(outd,"\n %4.0f",1800.0);
fprintf(outd,"\n %d",1000);
fprintf(outd,"\n %d",100);
fprintf(outd,"\n %d",100);
fprintf(outd,"\n %d",0);
fprintf(outd,"\n %d",10);
fprintf(outd,"\n %3.1f",0.0);
fprintf(outd,"\n %lf",6.0);
fprintf(outd,"\n %d",Atoms);
fprintf(outd,"\n %d",(int)size_x*size_y);
fprintf(outd,"\n %d",NUMBER_OF_FROZEN_LAYERS);
fprintf(outd,"\n %d",INIT_LAYERS-NUMBER_OF_FROZEN_LAYERS);
switch (surface)
{
case 0:
fprintf(outd,"\n [001]");
break;
case 1:
fprintf(outd,"\n [011]");
break;
case 2:
fprintf(outd,"\n [111]");
break;
default:
printf("Too Strange, too much different surface=%d \n", surface);
exit(0);
break;
}
fclose(outd);
}
comments(fn);
system("more -d aaa_info.txt");
exit(0);
}
}
//=============================//
void comments(char *fn)
{
FILE *outd;
if((outd=fopen(fn,"a"))!=NULL)
{
fprintf(outd,"\n *************** Main Parameters ***************** ");
fprintf(outd,"\n Regime: 0- start, 1 - equilibration, 2- measurements");
fprintf(outd,"\n Temperature");
fprintf(outd,"\n Number of steps in order to approach to the equilibrium ");
fprintf(outd,"\n Number of steps between two succesive measurements ");
fprintf(outd,"\n Number of measurements ");
fprintf(outd,"\n Initial step ");
fprintf(outd,"\n Number of the aviz pictures");
fprintf(outd,"\n Percent of the accomplished job");
fprintf(outd,"\n NVE/NVT parameter");
fprintf(outd,"\n Number of the atoms ");
fprintf(outd,"\n Number of the atoms in a layer ");
fprintf(outd,"\n Number of frozen layers");
fprintf(outd,"\n Number of free layers");
fprintf(outd,"\n The surface type: ");
fprintf(outd,"\n *************** Main Parameters ***************** ");
fclose(outd);
}
else
{ fprintf(stderr,"\n Can't open the file:\"%s\"",fn); }
}
//====================================//
void backup(char *fn,int mode)
{
extern long step;
if(mode)
{
save_configuration(mode) ;
save_inform(fn);
}
}
//====================================//
void add_info()
{
extern double lx,ly,lz,a0;
fprintf(stderr,"\n \n \t Lattice constants: a0(0)=%g a0(T)=%g",a00,a0);
fprintf(stderr,"\n \t Lattice box: [X Y Z] =%g x %g x %g ",lx,ly,lz);
}
//====================================//
void aviz_layer(char *fn)
{
extern double r[Atoms][3];
extern layers[NUMBER_OF_LAYERS][Atoms];
FILE * fd;
int j,i,k;
char *color[]=
{"a1","a2","a3","a4","a5",
"a6","a7","a8","a9","b1",
"b2","b3","b4","b5","b6",
"b7","b8","b9","c1","c2",
"c3","c4","c5","c6","c7",
"c8","c9","d1","d2","d3"} ;
if((fd=fopen(fn,"w"))!=NULL)
{ fprintf(fd," %d ", Atoms);
switch (surface)
{
case 0:
fprintf(fd,"\n ### Surface simulations (001) ");
break;
case 1:
fprintf(fd,"\n ### Surface simulations (011) ");
break;
default :
fprintf(fd,"\n ### Surface simulations (111) ");
break;
}
for (j=0;jFROZEN_ATOMS;j++)
fprintf(fd,"\n Cu %lf %lf %lf ", r[j][0],r[j][1],r[j][2]);
for(i=0;iNUMBER_OF_LAYERS;i++)
for (j=1;j=layers[i][0] ; j++)
{
k= layers[i][j];
fprintf(fd,"\n %s %lf %lf %lf ",color[i],r[k][0],r[k][1],r[k][2]);
}
fclose(fd);
}
else { fprintf(stderr,"I can't open the file %s !!!",fn); exit(0);}
}
*************** ``Makefile '' *****************
The Makefile for complitation and creation an executable file vanadium
INCS=initdata.h func_cp.h constants.h potential.h
CC=gcc -O3
LIBS=-lm
OBJS= measurements.o equilibration.o layers.o temp_prof.o bcc_setup.o compute_f.o density_profile.o evap_pto.o initialization.o md.o pair_correlation.o predictor_corrector.o read_save.o build_neigh_list.o order_parameter.o diffusion.o
vanadium: $(OBJS) $(INCS)
$(CC) $(OBJS) -o vanadium $(LIBS)
order_parameter.o: order_parameter.c $(INCS)
$(CC) -c order_parameter.c
pair_correlation.o: pair_correlation.c $(INCS)
$(CC) -c pair_correlation.c
bcc_setup.o: bcc_setup.c $(INCS)
$(CC) -c bcc_setup.c
build_neigh_list.o: build_neigh_list.c potential.h $(INCS)
$(CC) -c build_neigh_list.c
density_profile.o: density_profile.c $(INCS)
$(CC) -c density_profile.c
compute_f.o: compute_f.c $(INCS)
$(CC) -c compute_f.c
evap_pto.o: evap_pto.c $(INCS)
$(CC) -c evap_pto.c
initialization.o: initialization.c $(INCS)
$(CC) -c initialization.c
predictor_corrector.o: predictor_corrector.c $(INCS)
$(CC) -c predictor_corrector.c
read_save.o: read_save.c $(INCS)
$(CC) -c read_save.c
md.o: md.c $(INCS)
$(CC) -c md.c
diffusion.o: diffusion.c $(INCS)
$(CC) -c diffusion.c
temp_prof.o: temp_prof.c $(INCS)
$(CC) -c temp_prof.c
layers.o: layers.c $(INCS)
$(CC) -c layers.c
equilibration.o: equilibration.c $(INCS)
$(CC) -c equilibration.c
measurements.o: measurements.c $(INCS)
$(CC) -c measurements.c
//************** The m-script file ``setup_surf_geom.m''************//
a0=1.0;
%Type of the cell: 0- [001]; 1- [011]; 2- [111];
keys=1;
keys=input('Type of the cell: 0- [001]; 1- [011]; 2- [111]');
% Draw (1) or not to draw(0) the super cell points
draw_super_cell=0;
% Draw (1) or not to draw(0) the xy,xz,yz projections
draw_sections=0;
if(keys==0)
%%%%%%%%%%%%% [001] %%%%%%%%%%%%%%%%%%%%%%%%%%%%
A=eye(3); B=inv(A);
Nx=5; Ny=5; Nz=5;
Lx=a0*Nx; Ly=a0*Ny; Lz=a0*Nz;
%%%%%%%%%%%%% [001] %%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
if(keys==1)
%%%%%%%%%%%%% [011] %%%%%%%%%%%%%%%%%%%%%%%%%%%%
A=[-1./sqrt(2) 1./sqrt(2) 0
0 0 1.0
1./sqrt(2) 1./sqrt(2) 0];
Nx=5; Ny=5; Nz=1;
B=inv(A);
Lx=sqrt(2)*a0*Nx; Ly=a0*Ny; Lz=sqrt(1/2)*a0*Nz;
%%%%%%%%%%%%% [001] %%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
if(keys==2)
%%%%%%%%%%%%% [111] %%%%%%%%%%%%%%%%%%%%%%%%%%%%
A=[-1./sqrt(2.) 1./sqrt(2.) 0
-1./sqrt(6.) -1./sqrt(6.) sqrt(2./3.)
1./sqrt(3.) 1./sqrt(3.) 1./sqrt(3.0)];
B=inv(A);
Nx=6; Ny=6; Nz=18;
Nx=input('Input the number of Nx=Ny atoms ');
Ny=Nx;
Nz=input('Input the number of Nz atoms ');
Lx=(a0*sqrt(2))*Nx; Ly=(a0*sqrt(3/2))*Ny; Lz=a0/(2*sqrt(3))*Nz;
%%%%%%%%%%%%% [111] %%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
% Define dimension of the new surface:
Q=[0 0 0
Lx 0 0
0 Ly 0
0 0 Lz
Lx Ly 0
0 Ly Lz
Lx 0 Lz
Lx Ly Lz];
% Change the frame and find: the Limits
Xmax=-1000000;Xmin=1000000; Ymax=-1000000;Ymin=1000000; Zmax=-1000000;Zmin=1000000;
QQ=[];
for i=1:8
x=Q(i,1);y=Q(i,2);z=Q(i,3);
xx=B(1,1)*x+B(1,2)*y+B(1,3)*z;
yy=B(2,1)*x+B(2,2)*y+B(2,3)*z;
zz=B(3,1)*x+B(3,2)*y+B(3,3)*z;
if(xxXmax)Xmax=xx;end
if(yyYmax)Ymax=yy;end
if(zzZmax)Zmax=zz;end
if(xxXmin)Xmin=xx;end
if(yyYmin)Ymin=yy;end
if(zzZmin)Zmin=zz;end
QQ=[QQ;[xx,yy,zz]];
end
%Now, let's builda b.c.c crystal
ix=floor((Xmax-Xmin)/a0+0.5);
iy=floor((Ymax-Ymin)/a0+0.5);
iz=2*floor((Zmax-Zmin)/a0+0.5);
r=[];
pa=1;
for i=0:ix
for j=0:iy
for k=0:iz+2
r(pa,1)=Xmin+i*a0;
r(pa,2)=Ymin+j*a0;
r(pa,3)=Zmin+k*a0/2;
if(rem(k,2)==0)
r(pa,1)=r(pa,1)+a0/2;
r(pa,2)=r(pa,2)+a0/2;
end
pa=pa+1;
end
end
end
%Let's cut of the sample now
s=0;
R=[];RR=[];
epss=-0.00001;
for i=1:pa-1
x=r(i,1);y=r(i,2);z=r(i,3);
xx=A(1,1)*x+A(1,2)*y+A(1,3)*z;
yy=A(2,1)*x+A(2,2)*y+A(2,3)*z;
zz=A(3,1)*x+A(3,2)*y+A(3,3)*z;
%check the limits Lx,Ly,Lz
if(xx = epss & xx Lx+epss)
if(yy = epss & yy Ly+epss)
if( zz= epss & zz Lz+epss)
RR=[RR;[r(i,1),r(i,2),r(i,3)]];
R=[R;[xx,yy,zz]];
s=s+1;
end
end
end end
% Save the configurations:
fd=fopen('r00.log','w');
for i=1:s
fprintf(fd,'%4.8f %4.8f %4.8f \n',R(i,1),R(i,2),R(i,3));
end
fclose(fd);