small mods to support site energies and kinetic energies per atom

diff --git a/moldyn.c b/moldyn.c
index 433be68..7d27965 100644 (file)
--- 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;i<moldyn->count;i++)
-               moldyn->ekin+=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v));
+       for(i=0;i<moldyn->count;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 (file)
--- 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 (file)
--- 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 (file)
--- 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 (file)
--- 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 (file)
--- 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 (file)
--- 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 */