#include "../moldyn.h"
#include "../math/math.h"
-//#include "harmonic_oscillator.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_3dvec force,distance;
double d,f;
double sc,equi_dist;
+ double energy;
- params=moldyn->pot2b_params;
+ params=moldyn->pot_params;
sc=params->spring_constant;
equi_dist=params->equilibrium_distance;
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);
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;
+}
+