X-Git-Url: https://www.hackdaworld.org/gitweb/?a=blobdiff_plain;f=potentials%2Flennard_jones.c;h=2424b9851acf651eee69e7efa840b269c68086c4;hb=HEAD;hp=7d45a177134c13080cb5002ad474117f2450eab7;hpb=caa3bc828974c35df2462fde737c31c0a618ee4e;p=physik%2Fposic.git

diff --git a/potentials/lennard_jones.c b/potentials/lennard_jones.c
index 7d45a17..2424b98 100644
--- a/potentials/lennard_jones.c
+++ b/potentials/lennard_jones.c
@@ -16,8 +16,64 @@
 #include <math.h>
 
 #include "../moldyn.h"
-#inlcude "../math/math.h"
-//#include "lennard_jones.h"
+#include "../math/math.h"
+#include "lennard_jones.h"
+
+int lennard_jones_set_params(t_moldyn *moldyn,int element) {
+
+	t_lj_params *p;
+	t_atom a,b;
+
+	// set cutoff before parameters (actually only necessary for some pots)
+	if(moldyn->cutoff==0.0) {
+		printf("[lennard jones] WARNING: no cutoff!\n");
+		return -1;
+	}
+
+	/* atoms for correction term */
+	a.r.x=0.0; a.r.y=0.0; a.r.z=0.0;
+	b.r.x=0.0; b.r.y=0.0; b.r.z=moldyn->cutoff;
+
+	/* alloc mem for potential parameters */
+	moldyn->pot_params=malloc(sizeof(t_lj_params));
+	if(moldyn->pot_params==NULL) {
+		perror("[lennard jones] pot params alloc");
+		return -1;
+	}
+
+	/* these are now lennard jones parameters */
+	p=moldyn->pot_params;
+
+	switch(element) {
+		case SI:
+			/* type: silicon */
+			p->sigma6=LJ_SIGMA_SI*LJ_SIGMA_SI*LJ_SIGMA_SI;
+			p->sigma6+=p->sigma6;
+			p->sigma12=p->sigma6*p->sigma6;
+			p->epsilon4=4.0*LJ_EPSILON_SI;
+			p->uc=0.0; // calc it now!
+			lennard_jones(moldyn,&a,&b,0);
+			p->uc=moldyn->energy;
+			moldyn->energy=0.0;
+			break;
+		case C:
+			/* type: carbon */
+			p->sigma6=LJ_SIGMA_C*LJ_SIGMA_C*LJ_SIGMA_C;
+			p->sigma6+=p->sigma6;
+			p->sigma12=p->sigma6*p->sigma6;
+			p->epsilon4=4.0*LJ_EPSILON_C;
+			p->uc=0.0; // calc it now!
+			lennard_jones(moldyn,&a,&b,0);
+			p->uc=moldyn->energy;
+			moldyn->energy=0.0;
+			break;
+		default:
+			printf("[lennard jones] WARNING: element\n");
+			return -1;
+	}
+
+	return 0;
+}
 
 int lennard_jones(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 
@@ -25,8 +81,9 @@ 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->pot2b_params;
+	params=moldyn->pot_params;
 	eps=params->epsilon4;
 	sig6=params->sigma6;
 	sig12=params->sigma12;
@@ -35,12 +92,16 @@ int lennard_jones(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 
 	v3_sub(&distance,&(aj->r),&(ai->r));
 	if(bc) check_per_bound(moldyn,&distance);
-	d=v3_absolute_square(&distance);	/* 1/r^2 */
+	d=v3_absolute_square(&distance);	/* r^2 */
 	if(d<=moldyn->cutoff_square) {
 		d=1.0/d; 			/* 1/r^2 */
 		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;
@@ -58,3 +119,17 @@ int lennard_jones(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 	return 0;
 }
 
+int lennard_jones_check_2b_bond(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
+
+	 t_3dvec dist;
+	 double d;
+
+	v3_sub(&dist,&(aj->r),&(ai->r));
+	if(bc) check_per_bound(moldyn,&dist);
+	d=v3_absolute_square(&dist);
+
+	if(d>moldyn->cutoff_square)
+		return FALSE;
+
+	return TRUE;
+}