From: hackbard Date: Mon, 18 Dec 2006 09:30:22 +0000 (+0000) Subject: foo X-Git-Url: https://hackdaworld.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=2d562a1946322ae5b70bdce180321f32adbeab62;p=physik%2Fposic.git foo --- diff --git a/Makefile b/Makefile index 42ff2bd..39457f7 100644 --- a/Makefile +++ b/Makefile @@ -2,8 +2,9 @@ CC=gcc-3.4 #CC=gcc CFLAGS=-Wall CFLAGS+=-O3 -CFLAGS+=-g CFLAGS+=-ffloat-store +CFLAGS+=-g +CFLAGS+=-DDEBUG LDFLAGS=-lm OBJS=visual/visual.o random/random.o diff --git a/moldyn.c b/moldyn.c index db91859..faa6b7c 100644 --- a/moldyn.c +++ b/moldyn.c @@ -541,13 +541,13 @@ int scale_volume(t_moldyn *moldyn) { } /* 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 */ - 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; -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) @@ -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); + + /* 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"); @@ -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); + + /* 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"); @@ -1562,6 +1580,15 @@ if(ai==&(moldyn->atom[0])) { /* 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"); @@ -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 ^ */ + + /* 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"); diff --git a/sic.c b/sic.c index c4355b3..aabb79c 100644 --- 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_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 */ - 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, - &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, - &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 ! */ @@ -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); - //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"); @@ -161,13 +161,13 @@ int main(int argc,char **argv) { /* create the simulation schedule */ printf("[sic] adding schedule\n"); - moldyn_add_schedule(&md,10000,.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); - moldyn_set_log(&md,VISUAL_STEP,20); + moldyn_set_log(&md,VISUAL_STEP,1); /* * let's do the actual md algorithm now