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

[physik/posic.git] / potentials / lennard_jones.c

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;
+}