small mods to support site energies and kinetic energies per atom

[physik/posic.git] / moldyn.c

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;
 }