adapted all potential to new scheme + necessary mods to main code

diff --git a/moldyn.c b/moldyn.c
index def5079..a1bfcae 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -257,6 +257,14 @@ int set_potential(t_moldyn *moldyn,u8 type) {
                        moldyn->func3b_k2=albe_mult_3bp_k2;
                        moldyn->check_2b_bond=albe_mult_check_2b_bond;
                        break;
+               case MOLDYN_POTENTIAL_HO:
+                       moldyn->func2b=harmonic_oscillator;
+                       moldyn->check_2b_bond=harmonic_oscillator_check_2b_bond;
+                       break;
+               case MOLDYN_POTENTIAL_LJ:
+                       moldyn->func2b=lennard_jones;
+                       moldyn->check_2b_bond=lennard_jones_check_2b_bond;
+                       break;
                default:
                        printf("[moldyn] set potential: unknown type %02x\n",
                               type);

diff --git a/moldyn.h b/moldyn.h
index de6157c..18c09c1 100644 (file)
--- a/moldyn.h
+++ b/moldyn.h
@@ -52,9 +52,9 @@ typedef struct s_atom {
 #define ATOM_ATTR_VA   0x04    /* visualize this atom */
 #define ATOM_ATTR_VB   0x08    /* visualize the bond of this atom */
 
-#define ATOM_ATTR_1BP          0x10    /* single paricle potential */
-#define ATOM_ATTR_2BP          0x20    /* pair potential */
-#define ATOM_ATTR_3BP          0x40    /* 3 body potential */ 
+#define ATOM_ATTR_1BP  0x10    /* single paricle potential */
+#define ATOM_ATTR_2BP  0x20    /* pair potential */
+#define ATOM_ATTR_3BP  0x40    /* 3 body potential */ 
 
 /* cell lists */
 typedef struct s_linkcell {
@@ -314,12 +314,15 @@ typedef struct s_ba {
 #define SCALE_DIRECT                   'D'
 
 /*
- * potential related phsical values / constants
- *
+ * usefull constants
  */
 
 #define ONE_THIRD              (1.0/3.0)
 
+/*
+ * element specific defines
+ */
+
 #define C                      0x06
 #define LC_C                   3.567                           /* A */
 #define M_C                    12.011                          /* amu */
@@ -330,11 +333,6 @@ typedef struct s_ba {
 
 #define LC_3C_SIC              4.3596                          /* A */
 
-#define LJ_SIGMA_SI            ((0.25*sqrt(3.0)*LC_SI)/1.122462)       /* A */
-//#define LJ_SIGMA_SI          (LC_SI/1.122462)                        /* A */
-//#define LJ_SIGMA_SI          (0.5*sqrt(2.0)*LC_SI/1.122462)          /* A */
-#define LJ_EPSILON_SI          (2.1678*EV)                             /* NA */
-
 /*
  * lattice types
  */

diff --git a/potentials/albe.c b/potentials/albe.c
index ddfb3c6..c01ab15 100644 (file)
--- a/potentials/albe.c
+++ b/potentials/albe.c
@@ -24,6 +24,12 @@ int albe_mult_set_params(t_moldyn *moldyn,int element1,int element2) {
 
        t_albe_mult_params *p;
 
+       // set cutoff before parameters (actually only necessary for some pots)
+       if(moldyn->cutoff==0.0) {
+               printf("[albe] WARNING: no cutoff!\n");
+               return -1;
+       }
+
        /* alloc mem for potential parameters */
        moldyn->pot_params=malloc(sizeof(t_albe_mult_params));
        if(moldyn->pot_params==NULL) {

diff --git a/potentials/harmonic_oscillator.c b/potentials/harmonic_oscillator.c
index b203453..611f518 100644 (file)
--- a/potentials/harmonic_oscillator.c
+++ b/potentials/harmonic_oscillator.c
@@ -19,6 +19,46 @@
 #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;
@@ -55,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;
+}
+

diff --git a/potentials/harmonic_oscillator.h b/potentials/harmonic_oscillator.h
index c9c9891..13125d9 100644 (file)
--- a/potentials/harmonic_oscillator.h
+++ b/potentials/harmonic_oscillator.h
@@ -15,6 +15,19 @@ typedef struct s_ho_params {
 } t_ho_params;
 
 /* function prototype */
+int harmonic_oscillator_set_params(t_moldyn *moldyn,int element);
 int harmonic_oscillator(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc);
+int harmonic_oscillator_check_2b_bond(t_moldyn *moldyn,
+                                      t_atom *ai,t_atom *aj,u8 bc);
+
+/* harmonic oscillator potential parameter defines */
+
+// silicon
+#define HO_SC_SI               1
+#define HO_ED_SI               (0.25*sqrt(3.0)*LC_SI)
+
+// carbon
+#define HO_SC_C                        1
+#define HO_ED_C                        (0.25*sqrt(3.0)*LC_C)
 
 #endif

diff --git a/potentials/lennard_jones.c b/potentials/lennard_jones.c
index 4ff1b43..2424b98 100644 (file)
--- a/potentials/lennard_jones.c
+++ b/potentials/lennard_jones.c
@@ -19,6 +19,62 @@
 #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) {
 
        t_lj_params *params;
@@ -63,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;
+}

diff --git a/potentials/lennard_jones.h b/potentials/lennard_jones.h
index 617658b..8c1c36b 100644 (file)
--- a/potentials/lennard_jones.h
+++ b/potentials/lennard_jones.h
@@ -17,6 +17,18 @@ typedef struct s_lj_params {
 } t_lj_params;
 
 /* function prototype */
+int lennard_jones_set_params(t_moldyn *moldyn,int element);
 int lennard_jones(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc);
+int lennard_jones_check_2b_bond(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc);
+
+/* lennard jones potential parameters */
+
+// silicon
+#define LJ_SIGMA_SI            ((0.25*sqrt(3.0)*LC_SI)/1.122462)
+#define LJ_EPSILON_SI          (2.1678*EV) // TODO
+
+// carbob
+#define LJ_SIGMA_C             ((0.25*sqrt(3.0)*LC_C)/1.122462)
+#define LJ_EPSILON_C           (2.1678*EV) // TODO
 
 #endif

diff --git a/potentials/tersoff.c b/potentials/tersoff.c
index b83bfd6..d322b69 100644 (file)
--- a/potentials/tersoff.c
+++ b/potentials/tersoff.c
@@ -24,6 +24,12 @@ int tersoff_mult_set_params(t_moldyn *moldyn,int element1,int element2) {
 
        t_tersoff_mult_params *p;
 
+       // set cutoff before parameters (actually only necessary for some pots)
+       if(moldyn->cutoff==0.0) {
+               printf("[tersoff] WARNING: no cutoff!\n");
+               return -1;
+       }
+
        /* alloc mem for potential parameters */
         moldyn->pot_params=malloc(sizeof(t_tersoff_mult_params));
         if(moldyn->pot_params==NULL) {

diff --git a/sic.c b/sic.c
index 0e14bb1..a3c4f74 100644 (file)
--- a/sic.c
+++ b/sic.c
@@ -262,9 +262,11 @@ int main(int argc,char **argv) {
 
        /* choose potential */
 #ifdef ALBE
-       set_potential(&md,MOLDYN_POTENTIAL_AM);
+       if(set_potential(&md,MOLDYN_POTENTIAL_AM)<0)
+               return -1;
 #else
-       set_potential(&md,MOLDYN_POTENTIAL_TM);
+       if(set_potential(&md,MOLDYN_POTENTIAL_TM)<0)
+               return -1;
 #endif
 
        /* cutoff radius & bondlen */
@@ -282,15 +284,17 @@ int main(int argc,char **argv) {
         * potential parameters
         */
 
+#ifndef ALBE
        /*
          * tersoff mult potential parameters for SiC
         */
        tersoff_mult_set_params(&md,SI,C);
-
+#else
        /*
          * albe mult potential parameters for SiC
         */
        albe_mult_set_params(&md,SI,C);
+#endif
 
        /* set (initial) dimensions of simulation volume */
 #ifdef ALBE