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:
150c5f8
)
foo
author
hackbard
<hackbard>
Mon, 18 Dec 2006 09:30:22 +0000
(09:30 +0000)
committer
hackbard
<hackbard>
Mon, 18 Dec 2006 09:30:22 +0000
(09:30 +0000)
Makefile
patch
|
blob
|
history
moldyn.c
patch
|
blob
|
history
sic.c
patch
|
blob
|
history
diff --git
a/Makefile
b/Makefile
index
42ff2bd
..
39457f7
100644
(file)
--- a/
Makefile
+++ b/
Makefile
@@
-2,8
+2,9
@@
CC=gcc-3.4
#CC=gcc
CFLAGS=-Wall
CFLAGS+=-O3
#CC=gcc
CFLAGS=-Wall
CFLAGS+=-O3
-CFLAGS+=-g
CFLAGS+=-ffloat-store
CFLAGS+=-ffloat-store
+CFLAGS+=-g
+CFLAGS+=-DDEBUG
LDFLAGS=-lm
OBJS=visual/visual.o random/random.o
LDFLAGS=-lm
OBJS=visual/visual.o random/random.o
diff --git
a/moldyn.c
b/moldyn.c
index
db91859
..
faa6b7c
100644
(file)
--- a/
moldyn.c
+++ b/
moldyn.c
@@
-541,13
+541,13
@@
int scale_volume(t_moldyn *moldyn) {
}
/* just a guess so far ... */
}
/* just a guess so far ... */
- v=
sqrt(virial.xx*virial.xx+virial.yy*virial.yy+virial.zz+virial.zz)
;
+ v=
virial.xx+virial.yy+virial.zz
;
printf("%f\n",v);
/* get pressure from virial */
printf("%f\n",v);
/* get pressure from virial */
- moldyn->p=moldyn->count*K_BOLTZMANN*moldyn->t
-
ONE_THIRD*v;
+ moldyn->p=moldyn->count*K_BOLTZMANN*moldyn->t
+
ONE_THIRD*v;
moldyn->p/=moldyn->volume;
moldyn->p/=moldyn->volume;
-printf("%f
\n",moldyn->p/(ATM)
);
+printf("%f
| %f\n",moldyn->p/(ATM),moldyn->p_ref/ATM
);
/* scale factor */
if(moldyn->pt_scale&P_SCALE_BERENDSEN)
/* scale factor */
if(moldyn->pt_scale&P_SCALE_BERENDSEN)
@@
-1427,6
+1427,15
@@
int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
/* add forces of 2bp (ij, ji) contribution
* dVij = dVji and we sum up both: no 1/2) */
v3_add(&(ai->f),&(ai->f),&force);
/* add forces of 2bp (ij, ji) contribution
* dVij = dVji and we sum up both: no 1/2) */
v3_add(&(ai->f),&(ai->f),&force);
+
+ /* virial */
+ ai->virial.xx-=force.x*dist_ij.x;
+ ai->virial.yy-=force.y*dist_ij.y;
+ ai->virial.zz-=force.z*dist_ij.z;
+ ai->virial.xy-=force.x*dist_ij.y;
+ ai->virial.xz-=force.x*dist_ij.z;
+ ai->virial.yz-=force.y*dist_ij.z;
+
#ifdef DEBUG
if(ai==&(moldyn->atom[0])) {
printf("dVij, dVji (2bp) contrib:\n");
#ifdef DEBUG
if(ai==&(moldyn->atom[0])) {
printf("dVij, dVji (2bp) contrib:\n");
@@
-1526,6
+1535,15
@@
int tersoff_mult_post_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
/* add force */
v3_add(&(ai->f),&(ai->f),&force);
/* add force */
v3_add(&(ai->f),&(ai->f),&force);
+
+ /* virial */
+ ai->virial.xx-=force.x*dist_ij->x;
+ ai->virial.yy-=force.y*dist_ij->y;
+ ai->virial.zz-=force.z*dist_ij->z;
+ ai->virial.xy-=force.x*dist_ij->y;
+ ai->virial.xz-=force.x*dist_ij->z;
+ ai->virial.yz-=force.y*dist_ij->z;
+
#ifdef DEBUG
if(ai==&(moldyn->atom[0])) {
printf("dVij (3bp) contrib:\n");
#ifdef DEBUG
if(ai==&(moldyn->atom[0])) {
printf("dVij (3bp) contrib:\n");
@@
-1562,6
+1580,15
@@
if(ai==&(moldyn->atom[0])) {
/* add force */
v3_add(&(ai->f),&(ai->f),&force);
/* add force */
v3_add(&(ai->f),&(ai->f),&force);
+
+ /* virial - plus sign, as dist_ij = - dist_ji - (really??) */
+ ai->virial.xx+=force.x*dist_ij->x;
+ ai->virial.yy+=force.y*dist_ij->y;
+ ai->virial.zz+=force.z*dist_ij->z;
+ ai->virial.xy+=force.x*dist_ij->y;
+ ai->virial.xz+=force.x*dist_ij->z;
+ ai->virial.yz+=force.y*dist_ij->z;
+
#ifdef DEBUG
if(ai==&(moldyn->atom[0])) {
printf("dVji (3bp) contrib:\n");
#ifdef DEBUG
if(ai==&(moldyn->atom[0])) {
printf("dVji (3bp) contrib:\n");
@@
-1830,6
+1857,15
@@
int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
v3_scale(&temp2,&temp2,tmp*B*exp(-mu*d_jk)*f_c_jk*0.5);
v3_add(&(ai->f),&(ai->f),&temp2); /* -1 skipped in f_a calc ^ */
/* scaled with 0.5 ^ */
v3_scale(&temp2,&temp2,tmp*B*exp(-mu*d_jk)*f_c_jk*0.5);
v3_add(&(ai->f),&(ai->f),&temp2); /* -1 skipped in f_a calc ^ */
/* scaled with 0.5 ^ */
+
+ /* virial */
+ ai->virial.xx-=temp2.x*dist_jk.x;
+ ai->virial.yy-=temp2.y*dist_jk.y;
+ ai->virial.zz-=temp2.z*dist_jk.z;
+ ai->virial.xy-=temp2.x*dist_jk.y;
+ ai->virial.xz-=temp2.x*dist_jk.z;
+ ai->virial.yz-=temp2.y*dist_jk.z;
+
#ifdef DEBUG
if(ai==&(moldyn->atom[0])) {
printf("dVjk (3bp) contrib:\n");
#ifdef DEBUG
if(ai==&(moldyn->atom[0])) {
printf("dVjk (3bp) contrib:\n");
diff --git
a/sic.c
b/sic.c
index
c4355b3
..
aabb79c
100644
(file)
--- a/
sic.c
+++ b/
sic.c
@@
-112,26
+112,26
@@
int main(int argc,char **argv) {
/* create the lattice / place atoms */
printf("[sic] creating atoms\n");
/* create the lattice / place atoms */
printf("[sic] creating atoms\n");
-
//
create_lattice(&md,DIAMOND,LC_SI,SI,M_SI,
-
//
ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
-
//
0,5,5,5);
-
//
moldyn_bc_check(&md);
+ create_lattice(&md,DIAMOND,LC_SI,SI,M_SI,
+ ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
+ 0,5,5,5);
+ moldyn_bc_check(&md);
/* testing configuration */
/* testing configuration */
-
r.x=2.8/2;
v.x=0;
-
r.y=0;
v.y=0;
-
r.z=0;
v.z=0;
- add_atom(&md,SI,M_SI,0,
-
ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP,//
|ATOM_ATTR_HB,
+
//r.x=2.8/2;
v.x=0;
+
//r.y=0;
v.y=0;
+
//r.z=0;
v.z=0;
+
//
add_atom(&md,SI,M_SI,0,
+
// ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP
|ATOM_ATTR_HB,
// ATOM_ATTR_2BP,
// ATOM_ATTR_2BP,
- &r,&v);
-
r.x=-2.8/2;
v.x=0;
-
r.y=0;
v.y=0;
-
r.z=0;
v.z=0;
- add_atom(&md,SI,M_SI,0,
-
ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP,//
|ATOM_ATTR_HB,
+
//
&r,&v);
+
//r.x=-2.8/2;
v.x=0;
+
//r.y=0;
v.y=0;
+
//r.z=0;
v.z=0;
+
//
add_atom(&md,SI,M_SI,0,
+
// ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP
|ATOM_ATTR_HB,
// ATOM_ATTR_2BP,
// ATOM_ATTR_2BP,
- &r,&v);
+
//
&r,&v);
/* setting a nearest neighbour distance for the moldyn checks */
set_nn_dist(&md,0.25*sqrt(3.0)*LC_SI); /* diamond ! */
/* setting a nearest neighbour distance for the moldyn checks */
set_nn_dist(&md,0.25*sqrt(3.0)*LC_SI); /* diamond ! */
@@
-153,7
+153,7
@@
int main(int argc,char **argv) {
printf("[sic] set p/t scaling\n");
//set_pt_scale(&md,P_SCALE_BERENDSEN,100.0,
// T_SCALE_BERENDSEN,100.0);
printf("[sic] set p/t scaling\n");
//set_pt_scale(&md,P_SCALE_BERENDSEN,100.0,
// T_SCALE_BERENDSEN,100.0);
-
//
set_pt_scale(&md,0,0,T_SCALE_BERENDSEN,100.0);
+ set_pt_scale(&md,0,0,T_SCALE_BERENDSEN,100.0);
/* initial thermal fluctuations of particles (in equilibrium) */
printf("[sic] thermal init\n");
/* initial thermal fluctuations of particles (in equilibrium) */
printf("[sic] thermal init\n");
@@
-161,13
+161,13
@@
int main(int argc,char **argv) {
/* create the simulation schedule */
printf("[sic] adding schedule\n");
/* create the simulation schedule */
printf("[sic] adding schedule\n");
- moldyn_add_schedule(&md,100
00,.1
);
+ moldyn_add_schedule(&md,100
,1.0
);
/* activate logging */
printf("[sic] activate logging\n");
moldyn_set_log_dir(&md,argv[1]);
moldyn_set_log(&md,LOG_TOTAL_ENERGY,1);
/* activate logging */
printf("[sic] activate logging\n");
moldyn_set_log_dir(&md,argv[1]);
moldyn_set_log(&md,LOG_TOTAL_ENERGY,1);
- moldyn_set_log(&md,VISUAL_STEP,
20
);
+ moldyn_set_log(&md,VISUAL_STEP,
1
);
/*
* let's do the actual md algorithm now
/*
* let's do the actual md algorithm now