]> hackdaworld.org Git - physik/posic.git/commitdiff
introduced albe_orig (much faster!) + small change for c2, d2, c2/d2 ...
authorhackbard <hackbard@sage.physik.uni-augsburg.de>
Wed, 30 Jul 2008 15:12:21 +0000 (17:12 +0200)
committerhackbard <hackbard@sage.physik.uni-augsburg.de>
Wed, 30 Jul 2008 15:12:21 +0000 (17:12 +0200)
Makefile
mdrun.c
moldyn.c
moldyn.h
potentials/albe_orig.c
potentials/albe_orig.h

index 47670a825c9f5cf3af48675999688c6945ebc9e4..1b19d50171120e788d68e0af830d6d2f00eceb9c 100644 (file)
--- a/Makefile
+++ b/Makefile
@@ -11,8 +11,9 @@ CFLAGS += -DALBE
 
 #CFLAGS += -DSTATIC_LISTS
 
+#CFLAGS += -DNDEBUG
 #CFLAGS += -DDEBUG
-#CFLAGS += -DDSTART=19 -DDEND=40 -DDATOM=5832
+#CFLAGS += -DDSTART=-1 -DDEND=3 -DDATOM=0
 #CFLAGS += -DVDEBUG
 
 LDFLAGS = -lm
@@ -20,7 +21,7 @@ LDFLAGS = -lm
 
 DEPS = moldyn.o random/random.o list/list.o
 DEPS += potentials/lennard_jones.o potentials/harmonic_oscillator.o
-DEPS += potentials/tersoff.o potentials/albe.o
+DEPS += potentials/tersoff.o potentials/albe.o potentials/albe_orig.o
 
 ALL = mdrun sic fluctuation_calc postproc pair_correlation_calc diffusion_calc
 ALL += bond_analyze search_bonds visual_atoms display_atom_data
diff --git a/mdrun.c b/mdrun.c
index 29d3160acdc9375864ff56788a293704e8c3ccee..907388aafbcc45f6291895d88ec960c4acc0fa83 100644 (file)
--- a/mdrun.c
+++ b/mdrun.c
@@ -11,6 +11,7 @@
 #include "potentials/harmonic_oscillator.h"
 #include "potentials/lennard_jones.h"
 #include "potentials/albe.h"
+#include "potentials/albe_orig.h"
 #ifdef TERSOFF_ORIG
 #include "potentials/tersoff_orig.h"
 #else
@@ -196,6 +197,8 @@ int mdrun_parse_config(t_mdrun *mdrun) {
                if(!strncmp(word[0],"potential",9)) {
                        if(!strncmp(word[1],"albe",4))
                                mdrun->potential=MOLDYN_POTENTIAL_AM;
+                       if(!strncmp(word[1],"albe_o",6))
+                               mdrun->potential=MOLDYN_POTENTIAL_AO;
                        if(!strncmp(word[1],"tersoff",7))
                                mdrun->potential=MOLDYN_POTENTIAL_TM;
                        if(!strncmp(word[1],"ho",2))
@@ -905,6 +908,11 @@ int main(int argc,char **argv) {
                                             mdrun.element1,
                                             mdrun.element2);
                        break;
+               case MOLDYN_POTENTIAL_AO:
+                       albe_orig_mult_set_params(&moldyn,
+                                                 mdrun.element1,
+                                                 mdrun.element2);
+                       break;
                case MOLDYN_POTENTIAL_TM:
                        tersoff_mult_set_params(&moldyn,
                                                mdrun.element1,
index 1aff5a8df207c2b5d22077c70ccbc31c128afdd6..440a7075d4fe5556e5a679c3a1fc2f8a033c1077 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -24,6 +24,7 @@
 #include "potentials/harmonic_oscillator.h"
 #include "potentials/lennard_jones.h"
 #include "potentials/albe.h"
+#include "potentials/albe_orig.h"
 #ifdef TERSOFF_ORIG
 #include "potentials/tersoff_orig.h"
 #else
@@ -222,15 +223,13 @@ int set_potential(t_moldyn *moldyn,u8 type) {
                        moldyn->func_j1_k1=tersoff_mult_3bp_k2;
                        moldyn->check_2b_bond=tersoff_mult_check_2b_bond;
                        break;
-               /*
-               case MOLDYN_POTENTIAL_AM:
-                       moldyn->func3b_j1=albe_mult_3bp_j1;
-                       moldyn->func3b_k1=albe_mult_3bp_k1;
-                       moldyn->func3b_j2=albe_mult_3bp_j2;
-                       moldyn->func3b_k2=albe_mult_3bp_k2;
-                       moldyn->check_2b_bond=albe_mult_check_2b_bond;
+               case MOLDYN_POTENTIAL_AO:
+                       moldyn->func_j1=albe_orig_mult_3bp_j1;
+                       moldyn->func_j1_k0=albe_orig_mult_3bp_k1;
+                       moldyn->func_j1c=albe_orig_mult_3bp_j2;
+                       moldyn->func_j1_k1=albe_orig_mult_3bp_k2;
+                       moldyn->check_2b_bond=albe_orig_mult_check_2b_bond;
                        break;
-               */
                case MOLDYN_POTENTIAL_AM:
                        moldyn->func_i0=albe_mult_i0;
                        moldyn->func_j0=albe_mult_i0_j0;
index 066ad3642cfa9768809c5c17f47001ebaa3aca71..b2d6a3405ff13d3f66a1d60557510ef4397aaa1f 100644 (file)
--- a/moldyn.h
+++ b/moldyn.h
@@ -307,6 +307,7 @@ typedef struct s_vb {
 #define MOLDYN_POTENTIAL_LJ            0x01
 #define MOLDYN_POTENTIAL_TM            0x02
 #define MOLDYN_POTENTIAL_AM            0x03
+#define MOLDYN_POTENTIAL_AO            0x04
 
 #define LOG_TOTAL_ENERGY               0x01
 #define LOG_TOTAL_MOMENTUM             0x02
index c01ab15106d51834980fed6f9c0836c9336efcb1..46e160630cf160dfec4a912f27796d444e84d2b6 100644 (file)
@@ -1,5 +1,5 @@
 /*
- * albe.c - albe potential
+ * albe_orig.c - albe potential
  *
  * author: Frank Zirkelbach <frank.zirkelbach@physik.uni-augsburg.de>
  *
 
 #include "../moldyn.h"
 #include "../math/math.h"
-#include "albe.h"
+#include "albe_orig.h"
 
 /* create mixed terms from parameters and set them */
-int albe_mult_set_params(t_moldyn *moldyn,int element1,int element2) {
+int albe_orig_mult_set_params(t_moldyn *moldyn,int element1,int element2) {
 
-       t_albe_mult_params *p;
+       t_albe_orig_mult_params *p;
 
        // set cutoff before parameters (actually only necessary for some pots)
        if(moldyn->cutoff==0.0) {
-               printf("[albe] WARNING: no cutoff!\n");
+               printf("[albe orig] WARNING: no cutoff!\n");
                return -1;
        }
 
        /* alloc mem for potential parameters */
        moldyn->pot_params=malloc(sizeof(t_albe_mult_params));
        if(moldyn->pot_params==NULL) {
-               perror("[albe] pot params alloc");
+               perror("[albe orig] pot params alloc");
                return -1;
        }
 
@@ -53,7 +53,10 @@ int albe_mult_set_params(t_moldyn *moldyn,int element1,int element2) {
                        p->mu[0]=ALBE_MU_SI;
                        p->gamma[0]=ALBE_GAMMA_SI;
                        p->c[0]=ALBE_C_SI;
+                       p->c2[0]=p->c[0]*p->c[0];
                        p->d[0]=ALBE_D_SI;
+                       p->d2[0]=p->d[0]*p->d[0];
+                       p->c2d2[0]=p->c2[0]/p->d2[0];
                        p->h[0]=ALBE_H_SI;
                        switch(element2) {
                                case C:
@@ -67,7 +70,10 @@ int albe_mult_set_params(t_moldyn *moldyn,int element1,int element2) {
                                        p->mu[1]=ALBE_MU_C;
                                        p->gamma[1]=ALBE_GAMMA_C;
                                        p->c[1]=ALBE_C_C;
+                                       p->c2[1]=p->c[1]*p->c[1];
                                        p->d[1]=ALBE_D_C;
+                                       p->d2[1]=p->d[1]*p->d[1];
+                                       p->c2d2[1]=p->c2[1]/p->d2[1];
                                        p->h[1]=ALBE_H_C;
                                        /* mixed type: silicon carbide */
                                        p->Smixed=ALBE_S_SIC;
@@ -79,25 +85,29 @@ int albe_mult_set_params(t_moldyn *moldyn,int element1,int element2) {
                                        p->mu_m=ALBE_MU_SIC;
                                        p->gamma_m=ALBE_GAMMA_SIC;
                                        p->c_mixed=ALBE_C_SIC;
+                                       p->c2_m=p->c_mixed*p->c_mixed;
                                        p->d_mixed=ALBE_D_SIC;
+                                       p->d2_m=p->d_mixed*p->d_mixed;
+                                       p->c2d2_m=p->c2_m/p->d2_m;
                                        p->h_mixed=ALBE_H_SIC;
                                        break;
                                default:
-                                       printf("[albe] WARNING: element2\n");
+                                       printf("[albe orig] WARNING: element2");
+                                       printf("\n");
                                        return -1;
                        }
                        break;
                default:
-                       printf("[albe] WARNING: element1\n");
+                       printf("[albe orig] WARNING: element1\n");
                        return -1;
        }
 
-       printf("[albe] parameter completion\n");
+       printf("[albe orig] parameter completion\n");
        p->S2[0]=p->S[0]*p->S[0];
        p->S2[1]=p->S[1]*p->S[1];
        p->S2mixed=p->Smixed*p->Smixed;
 
-       printf("[albe] mult parameter info:\n");
+       printf("[albe orig] mult parameter info:\n");
        printf("  S (A)  | %f | %f | %f\n",p->S[0],p->S[1],p->Smixed);
        printf("  R (A)  | %f | %f | %f\n",p->R[0],p->R[1],p->Rmixed);
        printf("  A (eV) | %f | %f | %f\n",p->A[0]/EV,p->A[1]/EV,p->Amixed/EV);
@@ -105,19 +115,19 @@ int albe_mult_set_params(t_moldyn *moldyn,int element1,int element2) {
        printf("  lambda | %f | %f | %f\n",p->lambda[0],p->lambda[1],
                                          p->lambda_m);
        printf("  mu     | %f | %f | %f\n",p->mu[0],p->mu[1],p->mu_m);
-       printf("  gamma  | %f | %f\n",p->gamma[0],p->gamma[1]);
-       printf("  c      | %f | %f\n",p->c[0],p->c[1]);
-       printf("  d      | %f | %f\n",p->d[0],p->d[1]);
-       printf("  h      | %f | %f\n",p->h[0],p->h[1]);
+       printf("  gamma  | %f | %f | %f\n",p->gamma[0],p->gamma[1],p->gamma_m);
+       printf("  c      | %f | %f | %f\n",p->c[0],p->c[1],p->c_mixed);
+       printf("  d      | %f | %f | %f\n",p->d[0],p->d[1],p->d_mixed);
+       printf("  h      | %f | %f | %f\n",p->h[0],p->h[1],p->h_mixed);
 
        return 0;
 }
 
 /* albe 3 body potential function (first ij loop) */
-int albe_mult_3bp_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
+int albe_orig_mult_3bp_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 
-       t_albe_mult_params *params;
-       t_albe_exchange *exchange;
+       t_albe_orig_mult_params *params;
+       t_albe_orig_exchange *exchange;
        unsigned char brand;
        double S2;
        t_3dvec dist_ij;
@@ -167,11 +177,11 @@ int albe_mult_3bp_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 }
 
 /* albe 3 body potential function (first k loop) */
-int albe_mult_3bp_k1(t_moldyn *moldyn,
-                     t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
+int albe_orig_mult_3bp_k1(t_moldyn *moldyn,
+                          t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
 
-       t_albe_mult_params *params;
-       t_albe_exchange *exchange;
+       t_albe_orig_mult_params *params;
+       t_albe_orig_exchange *exchange;
        unsigned char brand;
        double R,S,S2;
        t_3dvec dist_ij,dist_ik;
@@ -180,11 +190,15 @@ int albe_mult_3bp_k1(t_moldyn *moldyn,
        double f_c_ik,df_c_ik;
        int kcount;
 
+       /* continue if aj equals ak */
+       if(aj==ak)
+               return 0;
+
        params=moldyn->pot_params;
        exchange=&(params->exchange);
        kcount=exchange->kcount;
 
-       if(kcount>ALBE_MAXN) {
+       if(kcount>ALBE_ORIG_MAXN) {
                printf("FATAL: neighbours = %d\n",kcount);
                printf("  -> %d %d %d\n",ai->tag,aj->tag,ak->tag);
        }
@@ -200,6 +214,9 @@ int albe_mult_3bp_k1(t_moldyn *moldyn,
                exchange->c_i=&(params->c[brand]);
                exchange->d_i=&(params->d[brand]);
                exchange->h_i=&(params->h[brand]);
+               exchange->ci2=&(params->c2[brand]);
+               exchange->di2=&(params->d2[brand]);
+               exchange->ci2di2=&(params->c2d2[brand]);
        }
        else {
                R=params->Rmixed;
@@ -210,10 +227,10 @@ int albe_mult_3bp_k1(t_moldyn *moldyn,
                exchange->c_i=&(params->c_mixed);
                exchange->d_i=&(params->d_mixed);
                exchange->h_i=&(params->h_mixed);
+               exchange->ci2=&(params->c2_m);
+               exchange->di2=&(params->d2_m);
+               exchange->ci2di2=&(params->c2d2_m);
        }
-       exchange->ci2=*(exchange->c_i)**(exchange->c_i);
-       exchange->di2=*(exchange->d_i)**(exchange->d_i);
-       exchange->ci2di2=exchange->ci2/exchange->di2;
 
        /* dist_ik, d_ik2 */
        v3_sub(&dist_ik,&(ak->r),&(ai->r));
@@ -242,9 +259,9 @@ int albe_mult_3bp_k1(t_moldyn *moldyn,
 
        /* g_ijk */
        h_cos=*(exchange->h_i)+cos_theta; // + in albe formalism
-       d2_h_cos2=exchange->di2+(h_cos*h_cos);
-       frac=exchange->ci2/d2_h_cos2;
-       g=*(exchange->gamma_i)*(1.0+exchange->ci2di2-frac);
+       d2_h_cos2=*(exchange->di2)+(h_cos*h_cos);
+       frac=*(exchange->ci2)/d2_h_cos2;
+       g=*(exchange->gamma_i)*(1.0+*(exchange->ci2di2)-frac);
        dg=2.0*frac**(exchange->gamma_i)*h_cos/d2_h_cos2; // + in albe f..
 
        /* zeta sum += f_c_ik * g_ijk */
@@ -275,10 +292,10 @@ int albe_mult_3bp_k1(t_moldyn *moldyn,
        return 0;
 }
 
-int albe_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
+int albe_orig_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 
-       t_albe_mult_params *params;
-       t_albe_exchange *exchange;
+       t_albe_orig_mult_params *params;
+       t_albe_orig_exchange *exchange;
        t_3dvec force;
        double f_a,df_a,b,db,f_c,df_c;
        double f_r,df_r;
@@ -388,11 +405,11 @@ if(moldyn->time>DSTART&&moldyn->time<DEND) {
 }
 
 /* albe 3 body potential function (second k loop) */
-int albe_mult_3bp_k2(t_moldyn *moldyn,
-                     t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
+int albe_orig_mult_3bp_k2(t_moldyn *moldyn,
+                          t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
 
-       t_albe_mult_params *params;
-       t_albe_exchange *exchange;
+       t_albe_orig_mult_params *params;
+       t_albe_orig_exchange *exchange;
        int kcount;
        t_3dvec dist_ik,dist_ij;
        double d_ik2,d_ik,d_ij2,d_ij;
@@ -405,6 +422,10 @@ int albe_mult_3bp_k2(t_moldyn *moldyn,
        t_3dvec dcosdrj,dcosdrk;
        t_3dvec force,tmp;
 
+       /* continue if aj equals ak */
+       if(aj==ak)
+               return 0;
+
        params=moldyn->pot_params;
        exchange=&(params->exchange);
        kcount=exchange->kcount;
@@ -521,9 +542,10 @@ if(moldyn->time>DSTART&&moldyn->time<DEND) {
 
 }
 
-int albe_mult_check_2b_bond(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,u8 bc) {
+int albe_orig_mult_check_2b_bond(t_moldyn *moldyn,
+                                 t_atom *itom,t_atom *jtom,u8 bc) {
 
-       t_albe_mult_params *params;
+       t_albe_orig_mult_params *params;
        t_3dvec dist;
        double d;
        u8 brand;
index 15386e76e5392f997e26e7e2639108cc24a985a6..a52c7b2656e45fd7cfded4973346a1b778da1196 100644 (file)
@@ -1,50 +1,53 @@
 /*
- * albe.h - albe potential header file
+ * albe_orig.h - albe potential header file
  *
  * author: Frank Zirkelbach <frank.zirkelbach@physik.uni-augsburg.de>
  *
  */
 
-#ifndef ALBE_H
-#define ALBE_H
+#ifndef ALBE_ORIG_H
+#define ALBE_ORIG_H
 
-#define ALBE_MAXN      16*27
+/* albe constants */
+#include "albe.h"
+
+#define ALBE_ORIG_MAXN 16*27
 
 /* albe exchange type */
-typedef struct s_albe_exchange {
+typedef struct s_albe_orig_exchange {
 
        t_3dvec dist_ij;
        double d_ij2;
        double d_ij;
 
-       t_3dvec dist_ik[ALBE_MAXN];
-       double d_ik2[ALBE_MAXN];
-       double d_ik[ALBE_MAXN];
+       t_3dvec dist_ik[ALBE_ORIG_MAXN];
+       double d_ik2[ALBE_ORIG_MAXN];
+       double d_ik[ALBE_ORIG_MAXN];
 
-       double f_c_ik[ALBE_MAXN];
-       double df_c_ik[ALBE_MAXN];
+       double f_c_ik[ALBE_ORIG_MAXN];
+       double df_c_ik[ALBE_ORIG_MAXN];
 
-       double g[ALBE_MAXN];
-       double dg[ALBE_MAXN];
-       double cos_theta[ALBE_MAXN];
+       double g[ALBE_ORIG_MAXN];
+       double dg[ALBE_ORIG_MAXN];
+       double cos_theta[ALBE_ORIG_MAXN];
 
        double *gamma_i;
        double *c_i;
        double *d_i;
        double *h_i;
 
-       double ci2;
-       double di2;
-       double ci2di2;
+       double *ci2;
+       double *di2;
+       double *ci2di2;
 
        double zeta_ij;
        double pre_dzeta;
 
        int kcount;
-} t_albe_exchange;
+} t_albe_orig_exchange;
 
 /* albe mult (2!) potential parameters */
-typedef struct s_albe_mult_params {
+typedef struct s_albe_orig_mult_params {
        double S[2];            /* albe cutoff radii */
        double S2[2];           /* albe cutoff radii squared */
        double R[2];            /* albe cutoff radii */
@@ -65,67 +68,33 @@ typedef struct s_albe_mult_params {
        double gamma[2];
        double gamma_m;
        double c[2];
+       double c2[2];
        double c_mixed;
+       double c2_m;
        double d[2];
+       double d2[2];
        double d_mixed;
+       double d2_m;
        double h[2];
        double h_mixed;
+       double c2d2[2];
+       double c2d2_m;
 
-       t_albe_exchange exchange;       /* exchange between 2bp and 3bp calc */
-} t_albe_mult_params;
+       t_albe_orig_exchange exchange;  /* exchange between 2bp and 3bp calc */
+} t_albe_orig_mult_params;
 
 /* function prototypes */
-int albe_mult_set_params(t_moldyn *moldyn,int element1,int elemnt2);
-int albe_mult_3bp_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc);
-int albe_mult_3bp_k1(t_moldyn *moldyn,
-                        t_atom *ai,t_atom *aj,t_atom *ak,u8 bc);
-int albe_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc);
-int albe_mult_3bp_k2(t_moldyn *moldyn,
-                        t_atom *ai,t_atom *aj,t_atom *ak,u8 bc);
-int albe_mult_check_2b_bond(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,u8 bc);
+int albe_orig_mult_set_params(t_moldyn *moldyn,int element1,int elemnt2);
+int albe_orig_mult_3bp_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc);
+int albe_orig_mult_3bp_k1(t_moldyn *moldyn,
+                          t_atom *ai,t_atom *aj,t_atom *ak,u8 bc);
+int albe_orig_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc);
+int albe_orig_mult_3bp_k2(t_moldyn *moldyn,
+                          t_atom *ai,t_atom *aj,t_atom *ak,u8 bc);
+int albe_orig_mult_check_2b_bond(t_moldyn *moldyn,
+                                 t_atom *itom,t_atom *jtom,u8 bc);
 
 /* albe potential parameter defines */
-
-// silicon
-#define ALBE_R_SI              (2.82-0.14)
-#define ALBE_S_SI              (2.82+0.14)
-#define ALBE_A_SI              (3.24*EV/0.842)
-#define ALBE_B_SI              (-1.842*3.24*EV/0.842)
-#define ALBE_R0_SI             2.232
-#define ALBE_LAMBDA_SI         (1.4761*sqrt(2.0*1.842))
-#define ALBE_MU_SI             (1.4761*sqrt(2.0/1.842))
-#define ALBE_GAMMA_SI          0.114354
-#define ALBE_C_SI              2.00494
-#define ALBE_D_SI              0.81472
-#define ALBE_H_SI              0.259
-#define ALBE_LC_SI             5.429
-
-// carbon
-#define ALBE_R_C               (2.00-0.15)
-#define ALBE_S_C               (2.00+0.15)
-#define ALBE_A_C               (6.00*EV/1.167)
-#define ALBE_B_C               (-2.167*6.00*EV/1.167)
-#define ALBE_R0_C              1.4276
-#define ALBE_LAMBDA_C          (2.0099*sqrt(2.0*2.167))
-#define ALBE_MU_C              (2.0099*sqrt(2.0/2.167))
-#define ALBE_GAMMA_C           0.11233
-#define ALBE_C_C               181.910
-#define ALBE_D_C               6.28433
-#define ALBE_H_C               0.5556
-#define ALBE_LC_C              3.566
-
-// mixed: silicon carbide
-#define ALBE_R_SIC             (2.40-0.20)
-#define ALBE_S_SIC             (2.40+0.20)
-#define ALBE_A_SIC             (4.36*EV/0.847)
-#define ALBE_B_SIC             (-1.847*4.36*EV/0.847)
-#define ALBE_R0_SIC            1.79
-#define ALBE_LAMBDA_SIC                (1.6991*sqrt(2.0*1.847))
-#define ALBE_MU_SIC            (1.6991*sqrt(2.0/1.847))
-#define ALBE_GAMMA_SIC         0.011877
-#define ALBE_C_SIC             273987
-#define ALBE_D_SIC             180.314
-#define ALBE_H_SIC             0.68
-#define ALBE_LC_SIC            4.359
+// -> albe.h
 
 #endif