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:
7493af8
)
safety checkin
author
hackbard
<hackbard>
Wed, 13 Dec 2006 00:37:10 +0000
(
00:37
+0000)
committer
hackbard
<hackbard>
Wed, 13 Dec 2006 00:37:10 +0000
(
00:37
+0000)
math/math.h
patch
|
blob
|
history
moldyn.c
patch
|
blob
|
history
sic.c
patch
|
blob
|
history
diff --git
a/math/math.h
b/math/math.h
index
8ff1bec
..
147f4c0
100644
(file)
--- a/
math/math.h
+++ b/
math/math.h
@@
-80,7
+80,12
@@
static inline int v3_set(t_3dvec *vec,double *ptr) {
static inline int v3_copy(t_3dvec *trg,t_3dvec *src) {
static inline int v3_copy(t_3dvec *trg,t_3dvec *src) {
- memcpy(trg,src,sizeof(t_3dvec));
+ //memcpy(trg,src,sizeof(t_3dvec));
+
+ /* faster this way? */
+ trg->x=src->x;
+ trg->y=src->y;
+ trg->z=src->z;
return 0;
}
return 0;
}
diff --git
a/moldyn.c
b/moldyn.c
index
44b53b6
..
73becde
100644
(file)
--- a/
moldyn.c
+++ b/
moldyn.c
@@
-802,14
+802,16
@@
int velocity_verlet(t_moldyn *moldyn) {
v3_add(&(atom[i].r),&(atom[i].r),&delta);
v3_scale(&delta,&(atom[i].f),0.5*tau_square/atom[i].mass);
v3_add(&(atom[i].r),&(atom[i].r),&delta);
v3_add(&(atom[i].r),&(atom[i].r),&delta);
v3_scale(&delta,&(atom[i].f),0.5*tau_square/atom[i].mass);
v3_add(&(atom[i].r),&(atom[i].r),&delta);
+//if(i==5) printf("v: %f %f %f\n",atom[i].r.x,(atom[i].r.x+moldyn->dim.x/2)/moldyn->lc.x,2*atom[i].r.x/moldyn->dim.x);
check_per_bound(moldyn,&(atom[i].r));
check_per_bound(moldyn,&(atom[i].r));
+//if(i==5) printf("n: %f %f %f\n",atom[i].r.x,(atom[i].r.x+moldyn->dim.x/2)/moldyn->lc.x,2*atom[i].r.x/moldyn->dim.x);
/* velocities */
v3_scale(&delta,&(atom[i].f),0.5*tau/atom[i].mass);
v3_add(&(atom[i].v),&(atom[i].v),&delta);
}
/* velocities */
v3_scale(&delta,&(atom[i].f),0.5*tau/atom[i].mass);
v3_add(&(atom[i].v),&(atom[i].v),&delta);
}
-moldyn_bc_check(moldyn);
+
//
moldyn_bc_check(moldyn);
/* neighbour list update */
link_cell_update(moldyn);
/* neighbour list update */
link_cell_update(moldyn);
@@
-841,7
+843,6
@@
int potential_force_calc(t_moldyn *moldyn) {
t_linkcell *lc;
t_list neighbour_i[27];
t_list neighbour_i2[27];
t_linkcell *lc;
t_list neighbour_i[27];
t_list neighbour_i2[27];
- //t_list neighbour_j[27];
t_list *this,*that;
u8 bc_ij,bc_ik;
int dnlc;
t_list *this,*that;
u8 bc_ij,bc_ik;
int dnlc;
@@
-988,7
+989,6
@@
int check_per_bound(t_moldyn *moldyn,t_3dvec *a) {
if(a->z>=z) a->z-=dim->z;
else if(-a->z>z) a->z+=dim->z;
}
if(a->z>=z) a->z-=dim->z;
else if(-a->z>z) a->z+=dim->z;
}
-printf("%f %f %f\n",a->x,x,a->x/x);
return 0;
}
return 0;
}
@@
-1236,7
+1236,8
@@
int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
s_r=S-R;
arg=M_PI*(d_ij-R)/s_r;
f_c=0.5+0.5*cos(arg);
s_r=S-R;
arg=M_PI*(d_ij-R)/s_r;
f_c=0.5+0.5*cos(arg);
- df_c=-0.5*sin(arg)*(M_PI/(s_r*d_ij));
+ //df_c=-0.5*sin(arg)*(M_PI/(s_r*d_ij)); /* MARK! */
+ df_c=0.5*sin(arg)*(M_PI/(s_r*d_ij));
/* two body contribution (ij, ji) */
v3_scale(&force,&dist_ij,-df_c*f_r-df_r*f_c);
/* tell 3bp that S > r > R */
/* two body contribution (ij, ji) */
v3_scale(&force,&dist_ij,-df_c*f_r-df_r*f_c);
/* tell 3bp that S > r > R */
@@
-1507,7
+1508,8
@@
int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
s_r=S-R;
arg=M_PI*(d_ik-R)/s_r;
f_c_ik=0.5+0.5*cos(arg);
s_r=S-R;
arg=M_PI*(d_ik-R)/s_r;
f_c_ik=0.5+0.5*cos(arg);
- df_c_ik=-0.5*sin(arg)*(M_PI/(s_r*d_ik));
+ //df_c_ik=-0.5*sin(arg)*(M_PI/(s_r*d_ik)); /* MARK */
+ df_c_ik=0.5*sin(arg)*(M_PI/(s_r*d_ik));
/* zeta_ij */
exchange->zeta_ij+=f_c_ik*g;
/* zeta_ij */
exchange->zeta_ij+=f_c_ik*g;
@@
-1557,13
+1559,13
@@
int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
c2d2=exchange->cj2dj2;
/* cosine of theta_jik by scalaproduct */
c2d2=exchange->cj2dj2;
/* cosine of theta_jik by scalaproduct */
- rr=
v3_scalar_product(&dist_ij,&dist_jk); /* times -1
*/
+ rr=
-v3_scalar_product(&dist_ij,&dist_jk); /* -1, as ij -> ji
*/
dd=d_ij*d_jk;
cos_theta=rr/dd;
/* d_costheta */
dd=d_ij*d_jk;
cos_theta=rr/dd;
/* d_costheta */
- d_costheta1=1.0/
(d_jk*d_ij)
;
- d_costheta2=cos_theta/(d_ij*d_ij);
/* in fact -cos(), but ^ */
+ d_costheta1=1.0/
dd
;
+ d_costheta2=cos_theta/(d_ij*d_ij);
/* some usefull values */
h_cos=(h-cos_theta);
/* some usefull values */
h_cos=(h-cos_theta);
@@
-1579,6
+1581,9
@@
int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
v3_add(&temp1,&temp1,&temp2);
v3_scale(&temp1,&temp1,-2.0*frac*h_cos/d2_h_cos2); /* dg */
v3_add(&temp1,&temp1,&temp2);
v3_scale(&temp1,&temp1,-2.0*frac*h_cos/d2_h_cos2); /* dg */
+ /* store dg in temp2 and use it for dVjk later */
+ v3_copy(&temp2,&temp1);
+
/* f_c_jk + {d,}zeta contribution (df_c_jk = 0) */
dzeta=&(exchange->dzeta_ji);
if(d_jk<R) {
/* f_c_jk + {d,}zeta contribution (df_c_jk = 0) */
dzeta=&(exchange->dzeta_ji);
if(d_jk<R) {
@@
-1600,7
+1605,7
@@
int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
/* zeta_ji */
exchange->zeta_ji+=f_c_jk*g;
/* zeta_ji */
exchange->zeta_ji+=f_c_jk*g;
- /* dzeta_
ij
*/
+ /* dzeta_
ji
*/
v3_scale(&temp1,&temp1,f_c_jk);
v3_add(dzeta,dzeta,&temp1);
}
v3_scale(&temp1,&temp1,f_c_jk);
v3_add(dzeta,dzeta,&temp1);
}
@@
-1608,19
+1613,19
@@
int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
/* dV_jk stuff | add force contribution on atom i immediately */
if(exchange->d_ij_between_rs) {
zeta=f_c*g;
/* dV_jk stuff | add force contribution on atom i immediately */
if(exchange->d_ij_between_rs) {
zeta=f_c*g;
- v3_scale(&temp1,&temp
1
,f_c);
- v3_scale(&temp2,&dist_ij,df_c);
- v3_add(&temp
1,&temp1,&temp2);
+ v3_scale(&temp1,&temp
2
,f_c);
+ v3_scale(&temp2,&dist_ij,df_c
*g
);
+ v3_add(&temp
2,&temp2,&temp1); /* -> dzeta_jk in temp2 */
}
else {
zeta=g;
}
else {
zeta=g;
- // dzeta_jk is simply dg, which is
temp1
+ // dzeta_jk is simply dg, which is
stored in temp2
}
/* betajnj * zeta_jk ^ nj-1 */
tmp=exchange->betajnj*pow(zeta,(n-1.0));
}
/* betajnj * zeta_jk ^ nj-1 */
tmp=exchange->betajnj*pow(zeta,(n-1.0));
- tmp=-chi/2.0*pow(
1+tmp*zeta,-1.0/(2.0*n)-1
)*tmp;
- v3_scale(&temp
1,&temp1
,tmp*B*exp(-mu*d_jk)*f_c_jk*0.5);
- v3_add(&(ai->f),&(ai->f),&temp
1
); /* -1 skipped in f_a calc ^ */
+ tmp=-chi/2.0*pow(
(1+tmp*zeta),(-1.0/(2.0*n)-1)
)*tmp;
+ v3_scale(&temp
2,&temp2
,tmp*B*exp(-mu*d_jk)*f_c_jk*0.5);
+ v3_add(&(ai->f),&(ai->f),&temp
2
); /* -1 skipped in f_a calc ^ */
/* scaled with 0.5 ^ */
}
/* scaled with 0.5 ^ */
}
@@
-1637,14
+1642,37
@@
int moldyn_bc_check(t_moldyn *moldyn) {
t_atom *atom;
t_3dvec *dim;
int i;
t_atom *atom;
t_3dvec *dim;
int i;
+double x;
+u8 byte;
+int j,k;
atom=moldyn->atom;
dim=&(moldyn->dim);
atom=moldyn->atom;
dim=&(moldyn->dim);
+x=dim->x/2;
for(i=0;i<moldyn->count;i++) {
for(i=0;i<moldyn->count;i++) {
- if(atom[i].r.x>=dim->x/2||-atom[i].r.x>dim->x/2)
+ if(atom[i].r.x>=dim->x/2||-atom[i].r.x>dim->x/2)
{
printf("FATAL: atom %d: x: %.20f (%.20f)\n",
i,atom[i].r.x,dim->x/2);
printf("FATAL: atom %d: x: %.20f (%.20f)\n",
i,atom[i].r.x,dim->x/2);
+ printf("diagnostic:\n");
+ for(j=0;j<8;j++) {
+ memcpy(&byte,(u8 *)(&(atom[i].r.x))+j,1);
+ for(k=0;k<8;k++)
+ printf("%d%c",
+ ((byte)&(1<<k))?1:0,
+ (k==7)?'\n':'|');
+ }
+ printf("---------------\n");
+ for(j=0;j<8;j++) {
+ memcpy(&byte,(u8 *)(&x)+j,1);
+ for(k=0;k<8;k++)
+ printf("%d%c",
+ ((byte)&(1<<k))?1:0,
+ (k==7)?'\n':'|');
+ }
+ if(atom[i].r.x==x) printf("hier gleich!\n");
+ else printf("hier NICHT gleich!\n");
+ }
if(atom[i].r.y>=dim->y/2||-atom[i].r.y>dim->y/2)
printf("FATAL: atom %d: y: %.20f (%.20f)\n",
i,atom[i].r.y,dim->y/2);
if(atom[i].r.y>=dim->y/2||-atom[i].r.y>dim->y/2)
printf("FATAL: atom %d: y: %.20f (%.20f)\n",
i,atom[i].r.y,dim->y/2);
diff --git
a/sic.c
b/sic.c
index
5edb699
..
bc10e31
100644
(file)
--- a/
sic.c
+++ b/
sic.c
@@
-132,7
+132,8
@@
int main(int argc,char **argv) {
/* set temperature */
printf("[sic] setting temperature\n");
//set_temperature(&md,273.0+450.0);
/* set temperature */
printf("[sic] setting temperature\n");
//set_temperature(&md,273.0+450.0);
- set_temperature(&md,0.0);
+ set_temperature(&md,1.0);
+ //set_temperature(&md,273.0);
/* set p/t scaling */
printf("[sic] set p/t scaling\n");
/* set p/t scaling */
printf("[sic] set p/t scaling\n");