potentials/harmonic_oscillator.c

   1 /*
   2  * harmonic_oscillator.c - harmonic oscillator potential
   3  *
   4  * author: Frank Zirkelbach <frank.zirkelbach@physik.uni-augsburg.de>
   5  *
   6  */
   7
   8 #define _GNU_SOURCE
   9 #include <stdio.h>
  10 #include <stdlib.h>
  11 #include <string.h>
  12 #include <sys/types.h>
  13 #include <sys/stat.h>
  14 #include <fcntl.h>
  15 #include <unistd.h>
  16 #include <math.h>
  17
  18 #include "../moldyn.h"
  19 #include "../math/math.h"
  20 #include "harmonic_oscillator.h"
  21
  22 int harmonic_oscillator_set_params(t_moldyn *moldyn,int element) {
  23
  24         t_ho_params *p;
  25
  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");
  29                 return -1;
  30         }
  31
  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");
  36                 return -1;
  37         }
  38
  39         /* these are now ho parameters */
  40         p=moldyn->pot_params;
  41
  42         switch(element) {
  43                 case SI:
  44                         /* type: silicon */
  45                         p->spring_constant=HO_SC_SI;
  46                         p->equilibrium_distance=HO_ED_SI;
  47                         break;
  48                 case C:
  49                         /* type: carbon */
  50                         p->spring_constant=HO_SC_C;
  51                         p->equilibrium_distance=HO_ED_C;
  52                         break;
  53                 default:
  54                         printf("[harmonic oscillator] WARNING: element\n");
  55                         return -1;
  56         }
  57
  58         return 0;
  59 }
  60
  61 int harmonic_oscillator(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
  62
  63         t_ho_params *params;
  64         t_3dvec force,distance;
  65         double d,f;
  66         double sc,equi_dist;
  67         double energy;
  68
  69         params=moldyn->pot_params;
  70         sc=params->spring_constant;
  71         equi_dist=params->equilibrium_distance;
  72
  73         if(ai<aj) return 0;
  74
  75         v3_sub(&distance,&(aj->r),&(ai->r));
  76
  77         if(bc) check_per_bound(moldyn,&distance);
  78         d=v3_norm(&distance);
  79         if(d<=moldyn->cutoff) {
  80                 energy=(0.5*sc*(d-equi_dist)*(d-equi_dist));
  81                 moldyn->energy+=energy;
  82                 ai->e+=0.5*energy;
  83                 aj->e+=0.5*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);
  92         }
  93
  94         return 0;
  95 }
  96
  97 int harmonic_oscillator_check_2b_bond(t_moldyn *moldyn,
  98                                       t_atom *ai,t_atom *aj,u8 bc) {
  99
 100         t_3dvec distance;
 101         double d;
 102
 103         v3_sub(&distance,&(aj->r),&(ai->r));
 104         if(bc) check_per_bound(moldyn,&distance);
 105         d=v3_norm(&distance);
 106
 107         if(d>moldyn->cutoff)
 108                 return FALSE;
 109
 110         return TRUE;
 111 }
 112