X-Git-Url: https://www.hackdaworld.org/gitweb/?a=blobdiff_plain;f=potentials%2Fharmonic_oscillator.c;h=611f518b3eeba9eeed5d9323685652b435ae10e2;hb=155e1cfea83209d09c2a06ae4fb7f5e1652fc00a;hp=e2a80e4f3401e3d7562fa5152b99e7537a515336;hpb=296a35b943e922173ce648ec76a4472e287af108;p=physik%2Fposic.git

diff --git a/potentials/harmonic_oscillator.c b/potentials/harmonic_oscillator.c
index e2a80e4..611f518 100644
--- a/potentials/harmonic_oscillator.c
+++ b/potentials/harmonic_oscillator.c
@@ -19,12 +19,53 @@
 #include "../math/math.h"
 #include "harmonic_oscillator.h"
 
+int harmonic_oscillator_set_params(t_moldyn *moldyn,int element) {
+
+	t_ho_params *p;
+
+	// set cutoff before parameters (actually only necessary for some pots)
+	if(moldyn->cutoff==0.0) {
+		printf("[harmonic oscillator] WARNING: no cutoff!\n");
+		return -1;
+	}
+
+	/* alloc mem for potential parameters */
+	moldyn->pot_params=malloc(sizeof(t_ho_params));
+	if(moldyn->pot_params==NULL) {
+		perror("[harmonic oscillator] pot params alloc");
+		return -1;
+	}
+
+	/* these are now ho parameters */
+	p=moldyn->pot_params;
+
+	switch(element) {
+		case SI:
+			/* type: silicon */
+			p->spring_constant=HO_SC_SI;
+			p->equilibrium_distance=HO_ED_SI;
+			break;
+		case C:
+			/* type: carbon */
+			p->spring_constant=HO_SC_C;
+			p->equilibrium_distance=HO_ED_C;
+			break;
+		default:
+			printf("[harmonic oscillator] WARNING: element\n");
+			return -1;
+	}
+
+
+	return 0;
+}
+
 int harmonic_oscillator(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 
 	t_ho_params *params;
 	t_3dvec force,distance;
 	double d,f;
 	double sc,equi_dist;
+	double energy;
 
 	params=moldyn->pot_params;
 	sc=params->spring_constant;
@@ -37,7 +78,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);
@@ -51,3 +95,19 @@ int harmonic_oscillator(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 	return 0;
 }
 
+int harmonic_oscillator_check_2b_bond(t_moldyn *moldyn,
+                                      t_atom *ai,t_atom *aj,u8 bc) {
+
+	t_3dvec distance;
+	double d;
+
+	v3_sub(&distance,&(aj->r),&(ai->r));
+	if(bc) check_per_bound(moldyn,&distance);
+	d=v3_norm(&distance);
+
+	if(d>moldyn->cutoff)
+		return FALSE;
+
+	return TRUE;
+}
+