From 20409ee0c545235c9246edde2d3cda5de0ddabde Mon Sep 17 00:00:00 2001 From: hackbard Date: Thu, 19 Jul 2007 18:04:08 +0200 Subject: [PATCH] small mods to support site energies and kinetic energies per atom --- moldyn.c | 73 +++++++++++++++++--------------- moldyn.h | 6 ++- potentials/albe.c | 5 ++- potentials/harmonic_oscillator.c | 6 ++- potentials/lennard_jones.c | 6 ++- potentials/tersoff.c | 11 ++++- sic.c | 4 +- 7 files changed, 68 insertions(+), 43 deletions(-) diff --git a/moldyn.c b/moldyn.c index 433be68..7d27965 100644 --- a/moldyn.c +++ b/moldyn.c @@ -510,6 +510,7 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass, atom[ret].brand=brand; atom[ret].tag=count+ret; check_per_bound(moldyn,&(atom[ret].r)); + atom[ret].r_0=atom[ret].r; } /* update total system mass */ @@ -518,6 +519,40 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass, return ret; } +int add_atom(t_moldyn *moldyn,int element,double mass,u8 brand,u8 attr, + t_3dvec *r,t_3dvec *v) { + + t_atom *atom; + void *ptr; + int count; + + atom=moldyn->atom; + count=(moldyn->count)++; + + ptr=realloc(atom,(count+1)*sizeof(t_atom)); + if(!ptr) { + perror("[moldyn] realloc (add atom)"); + return -1; + } + moldyn->atom=ptr; + + atom=moldyn->atom; + atom[count].r=*r; + atom[count].v=*v; + atom[count].element=element; + atom[count].mass=mass; + atom[count].brand=brand; + atom[count].tag=count; + atom[count].attr=attr; + check_per_bound(moldyn,&(atom[count].r)); + atom[count].r_0=atom[count].r; + + /* update total system mass */ + total_mass_calc(moldyn); + + return 0; +} + /* cubic init */ int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) { @@ -630,38 +665,6 @@ int diamond_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) { return count; } -int add_atom(t_moldyn *moldyn,int element,double mass,u8 brand,u8 attr, - t_3dvec *r,t_3dvec *v) { - - t_atom *atom; - void *ptr; - int count; - - atom=moldyn->atom; - count=(moldyn->count)++; - - ptr=realloc(atom,(count+1)*sizeof(t_atom)); - if(!ptr) { - perror("[moldyn] realloc (add atom)"); - return -1; - } - moldyn->atom=ptr; - - atom=moldyn->atom; - atom[count].r=*r; - atom[count].v=*v; - atom[count].element=element; - atom[count].mass=mass; - atom[count].brand=brand; - atom[count].tag=count; - atom[count].attr=attr; - - /* update total system mass */ - total_mass_calc(moldyn); - - return 0; -} - int destroy_atoms(t_moldyn *moldyn) { if(moldyn->atom) free(moldyn->atom); @@ -1097,8 +1100,10 @@ double e_kin_calc(t_moldyn *moldyn) { atom=moldyn->atom; moldyn->ekin=0.0; - for(i=0;icount;i++) - moldyn->ekin+=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v)); + for(i=0;icount;i++) { + atom[i].ekin=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v)); + moldyn->ekin+=atom[i].ekin; + } return moldyn->ekin; } diff --git a/moldyn.h b/moldyn.h index 51ad482..a595585 100644 --- a/moldyn.h +++ b/moldyn.h @@ -33,11 +33,13 @@ typedef struct s_virial { /* the atom of the md simulation */ typedef struct s_atom { + t_3dvec r_0; /* initial position */ t_3dvec r; /* position */ t_3dvec v; /* velocity */ t_3dvec f; /* force */ t_virial virial; /* virial */ double e; /* site energy */ + double ekin; /* kinetic energy */ int element; /* number of element in pse */ double mass; /* atom mass */ u8 brand; /* brand id */ @@ -410,11 +412,11 @@ int moldyn_log_shutdown(t_moldyn *moldyn); int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass, u8 attr,u8 brand,int a,int b,int c,t_3dvec *origin); +int add_atom(t_moldyn *moldyn,int element,double mass,u8 brand,u8 attr, + t_3dvec *r,t_3dvec *v); int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin); int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin); int diamond_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin); -int add_atom(t_moldyn *moldyn,int element,double mass,u8 brand,u8 attr, - t_3dvec *r,t_3dvec *v); int destroy_atoms(t_moldyn *moldyn); int thermal_init(t_moldyn *moldyn,u8 equi_init); diff --git a/potentials/albe.c b/potentials/albe.c index c0a5fe4..28008b0 100644 --- a/potentials/albe.c +++ b/potentials/albe.c @@ -226,6 +226,7 @@ int albe_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { double d_ij,r0; unsigned char brand; double S,R,s_r,arg; + double energy; params=moldyn->pot_params; exchange=&(params->exchange); @@ -309,7 +310,9 @@ int albe_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { exchange->pre_dzeta=-0.5*f_a*f_c*db; /* energy contribution */ - moldyn->energy+=0.5*f_c*(f_r+b*f_a); + energy=0.5*f_c*(f_r+b*f_a); + moldyn->energy+=energy; + ai->e+=energy; /* reset k counter for second k loop */ exchange->kcount=0; diff --git a/potentials/harmonic_oscillator.c b/potentials/harmonic_oscillator.c index e2a80e4..b203453 100644 --- a/potentials/harmonic_oscillator.c +++ b/potentials/harmonic_oscillator.c @@ -25,6 +25,7 @@ int harmonic_oscillator(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { t_3dvec force,distance; double d,f; double sc,equi_dist; + double energy; params=moldyn->pot_params; sc=params->spring_constant; @@ -37,7 +38,10 @@ int harmonic_oscillator(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { if(bc) check_per_bound(moldyn,&distance); d=v3_norm(&distance); if(d<=moldyn->cutoff) { - moldyn->energy+=(0.5*sc*(d-equi_dist)*(d-equi_dist)); + energy=(0.5*sc*(d-equi_dist)*(d-equi_dist)); + moldyn->energy+=energy; + ai->e+=0.5*energy; + aj->e+=0.5*energy; /* f = -grad E; grad r_ij = -1 1/r_ij distance */ f=sc*(1.0-equi_dist/d); v3_scale(&force,&distance,f); diff --git a/potentials/lennard_jones.c b/potentials/lennard_jones.c index 1277c46..4ff1b43 100644 --- a/potentials/lennard_jones.c +++ b/potentials/lennard_jones.c @@ -25,6 +25,7 @@ int lennard_jones(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { t_3dvec force,distance; double d,h1,h2; double eps,sig6,sig12; + double energy; params=moldyn->pot_params; eps=params->epsilon4; @@ -41,7 +42,10 @@ int lennard_jones(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { h2=d*d; /* 1/r^4 */ h2*=d; /* 1/r^6 */ h1=h2*h2; /* 1/r^12 */ - moldyn->energy+=(eps*(sig12*h1-sig6*h2)-params->uc); + energy=(eps*(sig12*h1-sig6*h2)-params->uc); + moldyn->energy+=energy; /* total energy */ + ai->e+=0.5*energy; /* site energy */ + aj->e+=0.5*energy; h2*=d; /* 1/r^8 */ h1*=d; /* 1/r^14 */ h2*=6*sig6; diff --git a/potentials/tersoff.c b/potentials/tersoff.c index 6bfcf30..953b839 100644 --- a/potentials/tersoff.c +++ b/potentials/tersoff.c @@ -93,6 +93,7 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { int brand; double s_r; double arg; + double energy; printf("WARNING! - tersoff_mult_2bp is obsolete.\n"); printf("WARNING! - repulsive part handled in 3bp/j2 routine.\n"); @@ -173,7 +174,10 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { virial_calc(ai,&force,&dist_ij); /* energy 2bp contribution */ - moldyn->energy+=f_r*f_c; + energy=f_r*f_c; + moldyn->energy+=energy; + ai->e+=0.5*energy; + aj->e+=0.5*energy; return 0; } @@ -348,6 +352,7 @@ int tersoff_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { unsigned char brand; double ni,tmp; double S,R,s_r,arg; + double energy; params=moldyn->pot_params; exchange=&(params->exchange); @@ -436,7 +441,9 @@ int tersoff_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { exchange->pre_dzeta=-0.5*f_a*f_c*db; /* energy contribution */ - moldyn->energy+=0.5*f_c*(b*f_a+f_r); + energy=0.5*f_c*(b*f_a+f_r); + moldyn->energy+=energy; + ai->e+=energy; /* reset k counter for second k loop */ exchange->kcount=0; diff --git a/sic.c b/sic.c index 6dc0027..5f2a8db 100644 --- a/sic.c +++ b/sic.c @@ -21,11 +21,11 @@ #include "potentials/tersoff.h" #endif -#define INJECT 1600 +#define INJECT 1 #define NR_ATOMS 1 #define R_C 2.0 #define T_C 10.0 -#define LCNT 20 +#define LCNT 2 typedef struct s_hp { int a_count; /* atom count */ -- 2.20.1