first sic tests

diff --git a/Makefile b/Makefile
index 70fcb8b..ebad7e3 100644 (file)
--- a/Makefile
+++ b/Makefile
@@ -6,12 +6,13 @@ CFLAGS+=-ffloat-store
 CFLAGS+=-g
 #CFLAGS+=-DDEBUG
 #CFLAGS+=-DVDEBUG
+#CFLAGS+=-DTERSOFF_ORIG
 LDFLAGS=-lm
 
 SOURCE=moldyn.c visual/visual.c random/random.c list/list.c
 
-POT_SRC=potentials/tersoff.c
 #POT_SRC=potentials/tersoff_orig.c
+POT_SRC=potentials/tersoff.c
 POT_SRC+= potentials/lennard_jones.c potentials/harmonic_oscillator.c
 
 all: sic

diff --git a/moldyn.c b/moldyn.c
index 705a28c..070ce69 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -421,11 +421,11 @@ int moldyn_log_shutdown(t_moldyn *moldyn) {
  */
 
 int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
-                   u8 attr,u8 brand,int a,int b,int c) {
+                   u8 attr,u8 brand,int a,int b,int c,t_3dvec *origin) {
 
        int new,count;
        int ret;
-       t_3dvec origin;
+       t_3dvec orig;
        void *ptr;
        t_atom *atom;
 
@@ -447,24 +447,33 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
        atom=&(moldyn->atom[count]);
 
        /* no atoms on the boundaries (only reason: it looks better!) */
-       origin.x=0.5*lc;
-       origin.y=0.5*lc;
-       origin.z=0.5*lc;
+       if(!origin) {
+               orig.x=0.5*lc;
+               orig.y=0.5*lc;
+               orig.z=0.5*lc;
+       }
+       else {
+               orig.x=origin->x;
+               orig.y=origin->y;
+               orig.z=origin->z;
+       }
 
        switch(type) {
                case CUBIC:
                        set_nn_dist(moldyn,lc);
-                       ret=cubic_init(a,b,c,lc,atom,&origin);
+                       ret=cubic_init(a,b,c,lc,atom,&orig);
                        break;
                case FCC:
-                       v3_scale(&origin,&origin,0.5);
+                       if(!origin)
+                               v3_scale(&orig,&orig,0.5);
                        set_nn_dist(moldyn,0.5*sqrt(2.0)*lc);
-                       ret=fcc_init(a,b,c,lc,atom,&origin);
+                       ret=fcc_init(a,b,c,lc,atom,&orig);
                        break;
                case DIAMOND:
-                       v3_scale(&origin,&origin,0.25);
+                       if(!origin)
+                               v3_scale(&orig,&orig,0.25);
                        set_nn_dist(moldyn,0.25*sqrt(3.0)*lc);
-                       ret=diamond_init(a,b,c,lc,atom,&origin);
+                       ret=diamond_init(a,b,c,lc,atom,&orig);
                        break;
                default:
                        printf("unknown lattice type (%02x)\n",type);
@@ -1252,6 +1261,9 @@ int moldyn_integrate(t_moldyn *moldyn) {
 
        /* calculate initial forces */
        potential_force_calc(moldyn);
+#ifdef DEBUG
+return 0;
+#endif
 
        /* some stupid checks before we actually start calculating bullshit */
        if(moldyn->cutoff>0.5*moldyn->dim.x)

diff --git a/moldyn.h b/moldyn.h
index 5b12048..3bc40c0 100644 (file)
--- a/moldyn.h
+++ b/moldyn.h
@@ -242,12 +242,15 @@ typedef struct s_moldyn {
 #define ONE_THIRD              (1.0/3.0)
 
 #define C                      0x06
+#define LC_C                   (0.3567e-9*METER)               /* A */
 #define M_C                    12.011                          /* amu */
 
 #define SI                     0x0e
 #define LC_SI                  (0.543105e-9*METER)             /* A */
 #define M_SI                   28.08553                        /* amu */
 
+#define LC_3C_SIC              (0.43596e-9*METER)              /* A */
+
 #define LJ_SIGMA_SI            ((0.25*sqrt(3.0)*LC_SI)/1.122462)       /* A */
 //#define LJ_SIGMA_SI          (LC_SI/1.122462)                        /* A */
 //#define LJ_SIGMA_SI          (0.5*sqrt(2.0)*LC_SI/1.122462)                  /* A */
@@ -279,6 +282,8 @@ typedef struct s_moldyn {
 
 #define TM_CHI_SIC             0.9776
 
+#define TM_LC_3C_SIC           (0.432e-9*METER)                /* A */
+
 /*
  * lattice constants
  */
@@ -324,7 +329,7 @@ int moldyn_set_log(t_moldyn *moldyn,u8 type,int timer);
 int moldyn_log_shutdown(t_moldyn *moldyn);
 
 int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
-                   u8 attr,u8 brand,int a,int b,int c);
+                   u8 attr,u8 brand,int a,int b,int c,t_3dvec *origin);
 int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin);
 int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin);
 int diamond_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin);

diff --git a/potentials/tersoff.c b/potentials/tersoff.c
index a4a87e3..da92ea5 100644 (file)
--- a/potentials/tersoff.c
+++ b/potentials/tersoff.c
@@ -418,7 +418,8 @@ int tersoff_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 #endif
 
        /* virial */
-       virial_calc(ai,&force,&(exchange->dist_ij));
+       if(aj<ai)
+               virial_calc(ai,&force,&(exchange->dist_ij));
 
        /* dzeta prefactor = - 0.5 f_c f_a db */
        exchange->pre_dzeta=-0.5*f_a*f_c*db;
@@ -544,7 +545,8 @@ int tersoff_mult_3bp_k2(t_moldyn *moldyn,
 
        /* virial */
        //v3_scale(&force,&force,-1.0);
-       virial_calc(ai,&force,&dist_ij);
+       if(aj<ai)
+               virial_calc(ai,&force,&dist_ij);
 
        /* derivative wrt k */
        v3_scale(&force,&dist_ik,-1.0*dfcg); // dri rik = - drk rik
@@ -565,7 +567,8 @@ int tersoff_mult_3bp_k2(t_moldyn *moldyn,
 
        /* virial */
        //v3_scale(&force,&force,-1.0);
-       virial_calc(ai,&force,&dist_ik);
+       if(aj<ai)
+               virial_calc(ai,&force,&dist_ik);
        
        /* increase k counter */
        exchange->kcount++;     

diff --git a/potentials/tersoff_orig.c b/potentials/tersoff_orig.c
index 69b67ac..431cfeb 100644 (file)
--- a/potentials/tersoff_orig.c
+++ b/potentials/tersoff_orig.c
@@ -629,7 +629,7 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
                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 */
+               /* store dg for use in dV_jk */
                v3_copy(&temp2,&temp1);
 
                /* f_c_jk + {d,}zeta contribution (df_c_jk = 0) */
@@ -670,8 +670,19 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
                        // dzeta_jk is simply dg, which is stored in temp2
                }
                /* betajnj * zeta_jk ^ nj-1 */
+printf("FATAL db_jk calc!\n");
+printf("(z1 + z2 + z3 ...)^n != z1^n + z2^n + z3^n + ...\n");
+printf("st00pid me => tersoff_orig is obsolete!\n");
                tmp=exchange->betajnj*pow(zeta,(n-1.0));
                tmp=-chi/2.0*pow((1+tmp*zeta),(-1.0/(2.0*n)-1))*tmp;
+#ifdef DEBUG
+       if((ai->tag==0)&(aj->tag==864)) { // &(ak->tag==23)) {
+       printf("\n\n");
+       printf("db: ni zeta = %f %f\n",
+              n,zeta);
+       printf("\n\n");
+       }
+#endif
                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 ^ */

diff --git a/sic.c b/sic.c
index 2221025..3c815ba 100644 (file)
--- a/sic.c
+++ b/sic.c
@@ -13,8 +13,12 @@
 /* potential */
 #include "potentials/harmonic_oscillator.h"
 #include "potentials/lennard_jones.h"
+
+#ifdef TERSOFF_ORIG
+#include "potentials/tersoff_orig.h"
+#else
 #include "potentials/tersoff.h"
-//#include "potentials/tersoff_orig.h"
+#endif
 
 #define INJECT         1
 #define NR_ATOMS       4
@@ -97,14 +101,16 @@ int main(int argc,char **argv) {
 
        /* choose potential */
        set_potential1b(&md,tersoff_mult_1bp);
-       //set_potential2b(&md,tersoff_mult_2bp);
-       //set_potential3b_j1(&md,tersoff_mult_2bp);
-       //set_potential3b_k1(&md,tersoff_mult_3bp);
-       //set_potential3b_j3(&md,tersoff_mult_post_2bp);
+#ifdef TERSOFF_ORIG
+       set_potential3b_j1(&md,tersoff_mult_2bp);
+       set_potential3b_k1(&md,tersoff_mult_3bp);
+       set_potential3b_j2(&md,tersoff_mult_post_2bp);
+#else
        set_potential3b_j1(&md,tersoff_mult_3bp_j1);
        set_potential3b_k1(&md,tersoff_mult_3bp_k1);
        set_potential3b_j2(&md,tersoff_mult_3bp_j2);
        set_potential3b_k2(&md,tersoff_mult_3bp_k2);
+#endif
        //set_potential2b(&md,lennard_jones);
        //set_potential2b(&md,harmonic_oscillator);
        set_potential_params(&md,&tp);
@@ -164,7 +170,8 @@ int main(int argc,char **argv) {
        tersoff_mult_complete_params(&tp);
 
        /* set (initial) dimensions of simulation volume */
-       set_dim(&md,6*LC_SI,6*LC_SI,6*LC_SI,TRUE);
+       //set_dim(&md,6*LC_SI,6*LC_SI,6*LC_SI,TRUE);
+       set_dim(&md,6*TM_LC_3C_SIC,6*TM_LC_3C_SIC,6*TM_LC_3C_SIC,TRUE);
 
        /* set periodic boundary conditions in all directions */
        set_pbc(&md,TRUE,TRUE,TRUE);
@@ -172,10 +179,21 @@ int main(int argc,char **argv) {
        /* create the lattice / place atoms */
        //create_lattice(&md,CUBIC,LC_SI,SI,M_SI,
        //create_lattice(&md,FCC,LC_SI,SI,M_SI,
-       create_lattice(&md,DIAMOND,LC_SI,SI,M_SI,
-                      ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
+       //create_lattice(&md,DIAMOND,LC_SI,SI,M_SI,
+       //create_lattice(&md,DIAMOND,LC_C,C,M_C,
+       //               ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
        //               ATOM_ATTR_2BP|ATOM_ATTR_HB,
-                      0,6,6,6);
+       //               1,6,6,6,NULL);
+
+       /* create centered zinc blende lattice */
+       r.x=0.5*0.25*TM_LC_3C_SIC; r.y=r.x; r.z=r.x;
+       create_lattice(&md,FCC,TM_LC_3C_SIC,SI,M_SI,
+                      ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
+                      0,6,6,6,&r);
+       r.x+=0.25*TM_LC_3C_SIC; r.y=r.x; r.z=r.x;
+       create_lattice(&md,FCC,TM_LC_3C_SIC,C,M_C,
+                      ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
+                      1,6,6,6,&r);
        moldyn_bc_check(&md);
 
        /* testing configuration */
@@ -225,7 +243,7 @@ int main(int argc,char **argv) {
 
        /* create the simulation schedule */
        /* initial configuration */
-       moldyn_add_schedule(&md,10000,1.0);
+       moldyn_add_schedule(&md,1000,1.0);
        /* adding atoms */
        //for(inject=0;inject<INJECT;inject++) {
        //      /* injecting atom and run with enabled t scaling */
@@ -243,8 +261,8 @@ int main(int argc,char **argv) {
        moldyn_set_log(&md,LOG_TOTAL_ENERGY,1);
        moldyn_set_log(&md,LOG_TEMPERATURE,1);
        moldyn_set_log(&md,LOG_PRESSURE,1);
-       moldyn_set_log(&md,VISUAL_STEP,100);
-       moldyn_set_log(&md,SAVE_STEP,100);
+       moldyn_set_log(&md,VISUAL_STEP,10);
+       moldyn_set_log(&md,SAVE_STEP,10);
        moldyn_set_log(&md,CREATE_REPORT,0);
 
        /*
@@ -253,6 +271,9 @@ int main(int argc,char **argv) {
         * integration of newtons equations
         */
        moldyn_integrate(&md);
+#ifdef DEBUG
+return 0;
+#endif
 
        /* close */
        moldyn_shutdown(&md);