2 * harmonic_oscillator.c - harmonic oscillator potential
4 * author: Frank Zirkelbach <frank.zirkelbach@physik.uni-augsburg.de>
12 #include <sys/types.h>
18 #include "../moldyn.h"
19 #include "../math/math.h"
20 #include "harmonic_oscillator.h"
22 int harmonic_oscillator_set_params(t_moldyn *moldyn,int element) {
26 // set cutoff before parameters (actually only necessary for some pots)
27 if(moldyn->cutoff==0.0) {
28 printf("[harmonic oscillator] WARNING: no cutoff!\n");
32 /* alloc mem for potential parameters */
33 moldyn->pot_params=malloc(sizeof(t_ho_params));
34 if(moldyn->pot_params==NULL) {
35 perror("[harmonic oscillator] pot params alloc");
39 /* these are now ho parameters */
45 p->spring_constant=HO_SC_SI;
46 p->equilibrium_distance=HO_ED_SI;
50 p->spring_constant=HO_SC_C;
51 p->equilibrium_distance=HO_ED_C;
54 printf("[harmonic oscillator] WARNING: element\n");
61 int harmonic_oscillator(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
64 t_3dvec force,distance;
69 params=moldyn->pot_params;
70 sc=params->spring_constant;
71 equi_dist=params->equilibrium_distance;
75 v3_sub(&distance,&(aj->r),&(ai->r));
77 if(bc) check_per_bound(moldyn,&distance);
79 if(d<=moldyn->cutoff) {
80 energy=(0.5*sc*(d-equi_dist)*(d-equi_dist));
81 moldyn->energy+=energy;
84 /* f = -grad E; grad r_ij = -1 1/r_ij distance */
85 f=sc*(1.0-equi_dist/d);
86 v3_scale(&force,&distance,f);
87 v3_add(&(ai->f),&(ai->f),&force);
88 virial_calc(ai,&force,&distance);
89 virial_calc(aj,&force,&distance); /* f and d signe switched */
90 v3_scale(&force,&distance,-f);
91 v3_add(&(aj->f),&(aj->f),&force);
97 int harmonic_oscillator_check_2b_bond(t_moldyn *moldyn,
98 t_atom *ai,t_atom *aj,u8 bc) {
103 v3_sub(&distance,&(aj->r),&(ai->r));
104 if(bc) check_per_bound(moldyn,&distance);
105 d=v3_norm(&distance);