projects
/
physik
/
posic.git
/ commitdiff
commit
grep
author
committer
pickaxe
?
search:
re
summary
|
shortlog
|
log
|
commit
| commitdiff |
tree
raw
|
patch
|
inline
| side by side (parent:
2d562a1
)
bugfix dg of dV_ji,jk
author
hackbard
<hackbard>
Thu, 28 Dec 2006 14:15:13 +0000
(14:15 +0000)
committer
hackbard
<hackbard>
Thu, 28 Dec 2006 14:15:13 +0000
(14:15 +0000)
moldyn.c
patch
|
blob
|
history
moldyn.h
patch
|
blob
|
history
diff --git
a/moldyn.c
b/moldyn.c
index
faa6b7c
..
c21f59f
100644
(file)
--- a/
moldyn.c
+++ b/
moldyn.c
@@
-1129,6
+1129,9
@@
int potential_force_calc(t_moldyn *moldyn) {
}
#ifdef DEBUG
printf("\n\n");
}
#ifdef DEBUG
printf("\n\n");
+#endif
+#ifdef VDEBUG
+printf("\n\n");
#endif
return 0;
#endif
return 0;
@@
-1376,6
+1379,7
@@
int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
/* save for use in 3bp */
exchange->d_ij=d_ij;
/* save for use in 3bp */
exchange->d_ij=d_ij;
+ exchange->d_ij2=d_ij2;
exchange->dist_ij=dist_ij;
/* more constants */
exchange->dist_ij=dist_ij;
/* more constants */
@@
-1443,6
+1447,14
@@
if(ai==&(moldyn->atom[0])) {
printf("%f | %f\n",force.y,ai->f.y);
printf("%f | %f\n",force.z,ai->f.z);
}
printf("%f | %f\n",force.y,ai->f.y);
printf("%f | %f\n",force.z,ai->f.z);
}
+#endif
+#ifdef VDEBUG
+if(ai==&(moldyn->atom[0])) {
+ printf("dVij, dVji (2bp) contrib:\n");
+ printf("%f | %f\n",force.x*dist_ij.x,ai->virial.xx);
+ printf("%f | %f\n",force.y*dist_ij.y,ai->virial.yy);
+ printf("%f | %f\n",force.z*dist_ij.z,ai->virial.zz);
+}
#endif
/* energy 2bp contribution (ij, ji) is 0.5 f_r f_c ... */
#endif
/* energy 2bp contribution (ij, ji) is 0.5 f_r f_c ... */
@@
-1551,6
+1563,14
@@
if(ai==&(moldyn->atom[0])) {
printf("%f | %f\n",force.y,ai->f.y);
printf("%f | %f\n",force.z,ai->f.z);
}
printf("%f | %f\n",force.y,ai->f.y);
printf("%f | %f\n",force.z,ai->f.z);
}
+#endif
+#ifdef VDEBUG
+if(ai==&(moldyn->atom[0])) {
+ printf("dVij (3bp) contrib:\n");
+ printf("%f | %f\n",force.x*dist_ij->x,ai->virial.xx);
+ printf("%f | %f\n",force.y*dist_ij->y,ai->virial.yy);
+ printf("%f | %f\n",force.z*dist_ij->z,ai->virial.zz);
+}
#endif
/* add energy of 3bp sum */
#endif
/* add energy of 3bp sum */
@@
-1596,6
+1616,14
@@
if(ai==&(moldyn->atom[0])) {
printf("%f | %f\n",force.y,ai->f.y);
printf("%f | %f\n",force.z,ai->f.z);
}
printf("%f | %f\n",force.y,ai->f.y);
printf("%f | %f\n",force.z,ai->f.z);
}
+#endif
+#ifdef VDEBUG
+if(ai==&(moldyn->atom[0])) {
+ printf("dVji (3bp) contrib:\n");
+ printf("%f | %f\n",force.x*dist_ij->x,ai->virial.xx);
+ printf("%f | %f\n",force.y*dist_ij->y,ai->virial.yy);
+ printf("%f | %f\n",force.z*dist_ij->z,ai->virial.zz);
+}
#endif
return 0;
#endif
return 0;
@@
-1610,9
+1638,9
@@
int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
t_3dvec dist_ij,dist_ik,dist_jk;
t_3dvec temp1,temp2;
t_3dvec *dzeta;
t_3dvec dist_ij,dist_ik,dist_jk;
t_3dvec temp1,temp2;
t_3dvec *dzeta;
- double R,S,s_r;
+ double R,S,
S2,
s_r;
double B,mu;
double B,mu;
- double d_ij,d_ik,d_jk;
+ double d_ij,d_ik,d_jk
,d_ij2,d_ik2,d_jk2
;
double rr,dd;
double f_c,df_c;
double f_c_ik,df_c_ik,arg;
double rr,dd;
double f_c,df_c;
double f_c_ik,df_c_ik,arg;
@@
-1658,6
+1686,7
@@
int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
/* dist_ij, d_ij - this is < S_ij ! */
dist_ij=exchange->dist_ij;
d_ij=exchange->d_ij;
/* dist_ij, d_ij - this is < S_ij ! */
dist_ij=exchange->dist_ij;
d_ij=exchange->d_ij;
+ d_ij2=exchange->d_ij2;
/* f_c_ij, df_c_ij (same for ji) */
f_c=exchange->f_c;
/* f_c_ij, df_c_ij (same for ji) */
f_c=exchange->f_c;
@@
-1672,21
+1701,26
@@
int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
/* dist_ik, d_ik */
v3_sub(&dist_ik,&(ak->r),&(ai->r));
if(bc) check_per_bound(moldyn,&dist_ik);
/* dist_ik, d_ik */
v3_sub(&dist_ik,&(ak->r),&(ai->r));
if(bc) check_per_bound(moldyn,&dist_ik);
- d_ik
=v3_norm
(&dist_ik);
+ d_ik
2=v3_absolute_square
(&dist_ik);
/* ik constants */
brand=ai->brand;
if(brand==ak->brand) {
R=params->R[brand];
S=params->S[brand];
/* ik constants */
brand=ai->brand;
if(brand==ak->brand) {
R=params->R[brand];
S=params->S[brand];
+ S2=params->S2[brand];
}
else {
R=params->Rmixed;
S=params->Smixed;
}
else {
R=params->Rmixed;
S=params->Smixed;
+ S2=params->S2mixed;
}
/* zeta_ij/dzeta_ij contribution only for d_ik < S */
}
/* zeta_ij/dzeta_ij contribution only for d_ik < S */
- if(d_ik<S) {
+ if(d_ik2<S2) {
+
+ /* now we need d_ik */
+ d_ik=sqrt(d_ik2);
/* get constants_i from exchange data */
n=*(exchange->n_i);
/* get constants_i from exchange data */
n=*(exchange->n_i);
@@
-1704,8
+1738,8
@@
int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
/* d_costheta */
tmp=1.0/dd;
/* d_costheta */
tmp=1.0/dd;
- d_costheta1=cos_theta/
(d_ij*d_ij)
-tmp;
- d_costheta2=cos_theta/
(d_ik*d_ik)
-tmp;
+ d_costheta1=cos_theta/
d_ij2
-tmp;
+ d_costheta2=cos_theta/
d_ik2
-tmp;
/* some usefull values */
h_cos=(h-cos_theta);
/* some usefull values */
h_cos=(h-cos_theta);
@@
-1757,13
+1791,14
@@
int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
/* dist_jk, d_jk */
v3_sub(&dist_jk,&(ak->r),&(aj->r));
if(bc) check_per_bound(moldyn,&dist_jk);
/* dist_jk, d_jk */
v3_sub(&dist_jk,&(ak->r),&(aj->r));
if(bc) check_per_bound(moldyn,&dist_jk);
- d_jk
=v3_norm
(&dist_jk);
+ d_jk
2=v3_absolute_square
(&dist_jk);
/* jk constants */
brand=aj->brand;
if(brand==ak->brand) {
R=params->R[brand];
S=params->S[brand];
/* jk constants */
brand=aj->brand;
if(brand==ak->brand) {
R=params->R[brand];
S=params->S[brand];
+ S2=params->S2[brand];
B=params->B[brand];
mu=params->mu[brand];
chi=1.0;
B=params->B[brand];
mu=params->mu[brand];
chi=1.0;
@@
-1771,13
+1806,17
@@
int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
else {
R=params->Rmixed;
S=params->Smixed;
else {
R=params->Rmixed;
S=params->Smixed;
+ S2=params->S2mixed;
B=params->Bmixed;
mu=params->mu_m;
chi=params->chi;
}
/* zeta_ji/dzeta_ji contribution only for d_jk < S_jk */
B=params->Bmixed;
mu=params->mu_m;
chi=params->chi;
}
/* zeta_ji/dzeta_ji contribution only for d_jk < S_jk */
- if(d_jk<S) {
+ if(d_jk2<S2) {
+
+ /* now we need d_ik */
+ d_jk=sqrt(d_jk2);
/* constants_j from exchange data */
n=*(exchange->n_j);
/* constants_j from exchange data */
n=*(exchange->n_j);
@@
-1795,7
+1834,7
@@
int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
/* d_costheta */
d_costheta1=1.0/dd;
/* d_costheta */
d_costheta1=1.0/dd;
- d_costheta2=cos_theta/
(d_ij*d_ij)
;
+ d_costheta2=cos_theta/
d_ij2
;
/* some usefull values */
h_cos=(h-cos_theta);
/* some usefull values */
h_cos=(h-cos_theta);
@@
-1805,10
+1844,11
@@
int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
/* g(cos_theta) */
g=1.0+c2d2-frac;
/* g(cos_theta) */
g=1.0+c2d2-frac;
- /* d_costheta_
ij
and dg(cos_theta) - needed in any case! */
+ /* d_costheta_
jik
and dg(cos_theta) - needed in any case! */
v3_scale(&temp1,&dist_jk,d_costheta1);
v3_scale(&temp2,&dist_ij,-d_costheta2); /* ji -> ij => -1 */
v3_scale(&temp1,&dist_jk,d_costheta1);
v3_scale(&temp2,&dist_ij,-d_costheta2); /* ji -> ij => -1 */
- v3_add(&temp1,&temp1,&temp2);
+ //v3_add(&temp1,&temp1,&temp2);
+ v3_sub(&temp1,&temp1,&temp2); /* there is a minus! */
v3_scale(&temp1,&temp1,-2.0*frac*h_cos/d2_h_cos2); /* dg */
/* store dg in temp2 and use it for dVjk later */
v3_scale(&temp1,&temp1,-2.0*frac*h_cos/d2_h_cos2); /* dg */
/* store dg in temp2 and use it for dVjk later */
@@
-1873,6
+1913,14
@@
if(ai==&(moldyn->atom[0])) {
printf("%f | %f\n",temp2.y,ai->f.y);
printf("%f | %f\n",temp2.z,ai->f.z);
}
printf("%f | %f\n",temp2.y,ai->f.y);
printf("%f | %f\n",temp2.z,ai->f.z);
}
+#endif
+#ifdef VDEBUG
+if(ai==&(moldyn->atom[0])) {
+ printf("dVjk (3bp) contrib:\n");
+ printf("%f | %f\n",temp2.x*dist_jk.x,ai->virial.xx);
+ printf("%f | %f\n",temp2.y*dist_jk.y,ai->virial.yy);
+ printf("%f | %f\n",temp2.z*dist_jk.z,ai->virial.zz);
+}
#endif
}
#endif
}
diff --git
a/moldyn.h
b/moldyn.h
index
dd94b0e
..
23b5100
100644
(file)
--- a/
moldyn.h
+++ b/
moldyn.h
@@
-195,6
+195,7
@@
typedef struct s_tersoff_exchange {
double f_a,df_a;
t_3dvec dist_ij;
double f_a,df_a;
t_3dvec dist_ij;
+ double d_ij2;
double d_ij;
double chi;
double d_ij;
double chi;