#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;
t_3dvec force,distance;
double d,h1,h2;
double eps,sig6,sig12;
+ double energy;
params=moldyn->pot_params;
eps=params->epsilon4;
h2=d*d; /* 1/r^4 */
h2*=d; /* 1/r^6 */
h1=h2*h2; /* 1/r^12 */
- moldyn->energy+=(eps*(sig12*h1-sig6*h2)-params->uc);
+ energy=(eps*(sig12*h1-sig6*h2)-params->uc);
+ moldyn->energy+=energy; /* total energy */
+ ai->e+=0.5*energy; /* site energy */
+ aj->e+=0.5*energy;
h2*=d; /* 1/r^8 */
h1*=d; /* 1/r^14 */
h2*=6*sig6;
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;
+}