projects
/
physik
/
posic.git
/ blobdiff
commit
grep
author
committer
pickaxe
?
search:
re
summary
|
shortlog
|
log
|
commit
|
commitdiff
|
tree
raw
|
inline
| side by side
float-store bugfix (unstatisfying!) + slightly new design of sic.c
[physik/posic.git]
/
potentials
/
lennard_jones.c
diff --git
a/potentials/lennard_jones.c
b/potentials/lennard_jones.c
index
7d45a17
..
4ff1b43
100644
(file)
--- a/
potentials/lennard_jones.c
+++ b/
potentials/lennard_jones.c
@@
-16,8
+16,8
@@
#include <math.h>
#include "../moldyn.h"
#include <math.h>
#include "../moldyn.h"
-#in
lc
ude "../math/math.h"
-
//
#include "lennard_jones.h"
+#in
cl
ude "../math/math.h"
+#include "lennard_jones.h"
int lennard_jones(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
int lennard_jones(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
@@
-25,8
+25,9
@@
int lennard_jones(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
t_3dvec force,distance;
double d,h1,h2;
double eps,sig6,sig12;
t_3dvec force,distance;
double d,h1,h2;
double eps,sig6,sig12;
+ double energy;
- params=moldyn->pot
2b
_params;
+ params=moldyn->pot_params;
eps=params->epsilon4;
sig6=params->sigma6;
sig12=params->sigma12;
eps=params->epsilon4;
sig6=params->sigma6;
sig12=params->sigma12;
@@
-35,12
+36,16
@@
int lennard_jones(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
v3_sub(&distance,&(aj->r),&(ai->r));
if(bc) check_per_bound(moldyn,&distance);
v3_sub(&distance,&(aj->r),&(ai->r));
if(bc) check_per_bound(moldyn,&distance);
- d=v3_absolute_square(&distance); /*
1/
r^2 */
+ d=v3_absolute_square(&distance); /* r^2 */
if(d<=moldyn->cutoff_square) {
d=1.0/d; /* 1/r^2 */
h2=d*d; /* 1/r^4 */
if(d<=moldyn->cutoff_square) {
d=1.0/d; /* 1/r^2 */
h2=d*d; /* 1/r^4 */
+ h2*=d; /* 1/r^6 */
h1=h2*h2; /* 1/r^12 */
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;
h2*=d; /* 1/r^8 */
h1*=d; /* 1/r^14 */
h2*=6*sig6;