From: hackbard <hackbard@sage.physik.uni-augsburg.de>
Date: Thu, 19 Jul 2007 16:04:08 +0000 (+0200)
Subject: small mods to support site energies and kinetic energies per atom
X-Git-Url: https://hackdaworld.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=20409ee0c545235c9246edde2d3cda5de0ddabde;p=physik%2Fposic.git

small mods to support site energies and kinetic energies per atom
---

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;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
--- 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 */