X-Git-Url: https://hackdaworld.org/cgi-bin/gitweb.cgi?a=blobdiff_plain;f=potentials%2Ftersoff.c;h=16f771c7be4b614f10c6fd70b9476fb0bfd5e23f;hb=ea612b88a0588b8f46fafaebf3b37fb46c83c0cf;hp=87bb1874077c39ca7992030626de689ab1b0e57d;hpb=caa3bc828974c35df2462fde737c31c0a618ee4e;p=physik%2Fposic.git

diff --git a/potentials/tersoff.c b/potentials/tersoff.c
index 87bb187..16f771c 100644
--- a/potentials/tersoff.c
+++ b/potentials/tersoff.c
@@ -17,7 +17,7 @@
 
 #include "../moldyn.h"
 #include "../math/math.h"
-//#include "tersoff.h"
+#include "tersoff.h"
 
 /* create mixed terms from parameters and set them */
 int tersoff_mult_complete_params(t_tersoff_mult_params *p) {
@@ -59,7 +59,7 @@ int tersoff_mult_1bp(t_moldyn *moldyn,t_atom *ai) {
 	t_tersoff_exchange *exchange;
 	
 	brand=ai->brand;
-	params=moldyn->pot1b_params;
+	params=moldyn->pot_params;
 	exchange=&(params->exchange);
 
 	/*
@@ -95,7 +95,10 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 	double s_r;
 	double arg;
 
-	params=moldyn->pot2b_params;
+	/* use newtons third law */
+	//if(ai<aj) return 0;
+
+	params=moldyn->pot_params;
 	brand=aj->brand;
 	exchange=&(params->exchange);
 
@@ -119,27 +122,11 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 	 *
 	 */
 
-	/* constants */
-	if(brand==ai->brand) {
-		S=params->S[brand];
+	/* determine cutoff square */
+	if(brand==ai->brand)
 		S2=params->S2[brand];
-		R=params->R[brand];
-		A=params->A[brand];
-		B=params->B[brand];
-		lambda=params->lambda[brand];
-		mu=params->mu[brand];
-		exchange->chi=1.0;
-	}
-	else {
-		S=params->Smixed;
+	else
 		S2=params->S2mixed;
-		R=params->Rmixed;
-		A=params->Amixed;
-		B=params->Bmixed;
-		lambda=params->lambda_m;
-		mu=params->mu_m;
-		params->exchange.chi=params->chi;
-	}
 
 	/* dist_ij, d_ij */
 	v3_sub(&dist_ij,&(aj->r),&(ai->r));
@@ -166,12 +153,26 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 	exchange->d_j=&(params->d[brand]);
 	exchange->h_j=&(params->h[brand]);
 	if(brand==ai->brand) {
+		S=params->S[brand];
+		R=params->R[brand];
+		A=params->A[brand];
+		B=params->B[brand];
+		lambda=params->lambda[brand];
+		mu=params->mu[brand];
+		exchange->chi=1.0;
 		exchange->betajnj=exchange->betaini;
 		exchange->cj2=exchange->ci2;
 		exchange->dj2=exchange->di2;
 		exchange->cj2dj2=exchange->ci2di2;
 	}
 	else {
+		S=params->Smixed;
+		R=params->Rmixed;
+		A=params->Amixed;
+		B=params->Bmixed;
+		lambda=params->lambda_m;
+		mu=params->mu_m;
+		params->exchange.chi=params->chi;
 		exchange->betajnj=pow(*(exchange->beta_j),*(exchange->n_j));
 		exchange->cj2=params->c[brand]*params->c[brand];
 		exchange->dj2=params->d[brand]*params->d[brand];
@@ -210,12 +211,13 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 	v3_add(&(ai->f),&(ai->f),&force);
 
 	/* virial */
-	ai->virial.xx-=force.x*dist_ij.x;
-	ai->virial.yy-=force.y*dist_ij.y;
-	ai->virial.zz-=force.z*dist_ij.z;
-	ai->virial.xy-=force.x*dist_ij.y;
-	ai->virial.xz-=force.x*dist_ij.z;
-	ai->virial.yz-=force.y*dist_ij.z;
+	virial_calc(ai,&force,&dist_ij);
+	//ai->virial.xx-=force.x*dist_ij.x;
+	//ai->virial.yy-=force.y*dist_ij.y;
+	//ai->virial.zz-=force.z*dist_ij.z;
+	//ai->virial.xy-=force.x*dist_ij.y;
+	//ai->virial.xz-=force.x*dist_ij.z;
+	//ai->virial.yz-=force.y*dist_ij.z;
 
 #ifdef DEBUG
 if(ai==&(moldyn->atom[0])) {
@@ -281,7 +283,7 @@ int tersoff_mult_post_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 	double chi,ni,betaini,nj,betajnj;
 	double zeta;
 
-	params=moldyn->pot2b_params;
+	params=moldyn->pot_params;
 	exchange=&(params->exchange);
 
 	/* we do not run if f_c_ij was detected to be 0! */
@@ -325,12 +327,13 @@ int tersoff_mult_post_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 	v3_add(&(ai->f),&(ai->f),&force);
 
 	/* virial */
-	ai->virial.xx-=force.x*dist_ij->x;
-	ai->virial.yy-=force.y*dist_ij->y;
-	ai->virial.zz-=force.z*dist_ij->z;
-	ai->virial.xy-=force.x*dist_ij->y;
-	ai->virial.xz-=force.x*dist_ij->z;
-	ai->virial.yz-=force.y*dist_ij->z;
+	virial_calc(ai,&force,dist_ij);
+	//ai->virial.xx-=force.x*dist_ij->x;
+	//ai->virial.yy-=force.y*dist_ij->y;
+	//ai->virial.zz-=force.z*dist_ij->z;
+	//ai->virial.xy-=force.x*dist_ij->y;
+	//ai->virial.xz-=force.x*dist_ij->z;
+	//ai->virial.yz-=force.y*dist_ij->z;
 
 #ifdef DEBUG
 if(ai==&(moldyn->atom[0])) {
@@ -379,12 +382,13 @@ if(ai==&(moldyn->atom[0])) {
 
 	/* virial - plus sign, as dist_ij = - dist_ji - (really??) */
 // TEST ... with a minus instead
-	ai->virial.xx-=force.x*dist_ij->x;
-	ai->virial.yy-=force.y*dist_ij->y;
-	ai->virial.zz-=force.z*dist_ij->z;
-	ai->virial.xy-=force.x*dist_ij->y;
-	ai->virial.xz-=force.x*dist_ij->z;
-	ai->virial.yz-=force.y*dist_ij->z;
+	virial_calc(ai,&force,dist_ij);
+	//ai->virial.xx-=force.x*dist_ij->x;
+	//ai->virial.yy-=force.y*dist_ij->y;
+	//ai->virial.zz-=force.z*dist_ij->z;
+	//ai->virial.xy-=force.x*dist_ij->y;
+	//ai->virial.xz-=force.x*dist_ij->z;
+	//ai->virial.yz-=force.y*dist_ij->z;
 
 #ifdef DEBUG
 if(ai==&(moldyn->atom[0])) {
@@ -430,7 +434,7 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
 	double tmp;
 	int brand;
 
-	params=moldyn->pot3b_params;
+	params=moldyn->pot_params;
 	exchange=&(params->exchange);
 
 	if(!(exchange->run3bp))