From: hackbard Date: Wed, 16 Apr 2008 11:15:39 +0000 (+0200) Subject: adapted all potential to new scheme + necessary mods to main code X-Git-Url: https://hackdaworld.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=155e1cfea83209d09c2a06ae4fb7f5e1652fc00a;p=physik%2Fposic.git adapted all potential to new scheme + necessary mods to main code --- diff --git a/moldyn.c b/moldyn.c index def5079..a1bfcae 100644 --- a/moldyn.c +++ b/moldyn.c @@ -257,6 +257,14 @@ int set_potential(t_moldyn *moldyn,u8 type) { moldyn->func3b_k2=albe_mult_3bp_k2; moldyn->check_2b_bond=albe_mult_check_2b_bond; break; + case MOLDYN_POTENTIAL_HO: + moldyn->func2b=harmonic_oscillator; + moldyn->check_2b_bond=harmonic_oscillator_check_2b_bond; + break; + case MOLDYN_POTENTIAL_LJ: + moldyn->func2b=lennard_jones; + moldyn->check_2b_bond=lennard_jones_check_2b_bond; + break; default: printf("[moldyn] set potential: unknown type %02x\n", type); diff --git a/moldyn.h b/moldyn.h index de6157c..18c09c1 100644 --- a/moldyn.h +++ b/moldyn.h @@ -52,9 +52,9 @@ typedef struct s_atom { #define ATOM_ATTR_VA 0x04 /* visualize this atom */ #define ATOM_ATTR_VB 0x08 /* visualize the bond of this atom */ -#define ATOM_ATTR_1BP 0x10 /* single paricle potential */ -#define ATOM_ATTR_2BP 0x20 /* pair potential */ -#define ATOM_ATTR_3BP 0x40 /* 3 body potential */ +#define ATOM_ATTR_1BP 0x10 /* single paricle potential */ +#define ATOM_ATTR_2BP 0x20 /* pair potential */ +#define ATOM_ATTR_3BP 0x40 /* 3 body potential */ /* cell lists */ typedef struct s_linkcell { @@ -314,12 +314,15 @@ typedef struct s_ba { #define SCALE_DIRECT 'D' /* - * potential related phsical values / constants - * + * usefull constants */ #define ONE_THIRD (1.0/3.0) +/* + * element specific defines + */ + #define C 0x06 #define LC_C 3.567 /* A */ #define M_C 12.011 /* amu */ @@ -330,11 +333,6 @@ typedef struct s_ba { #define LC_3C_SIC 4.3596 /* 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 */ -#define LJ_EPSILON_SI (2.1678*EV) /* NA */ - /* * lattice types */ diff --git a/potentials/albe.c b/potentials/albe.c index ddfb3c6..c01ab15 100644 --- a/potentials/albe.c +++ b/potentials/albe.c @@ -24,6 +24,12 @@ int albe_mult_set_params(t_moldyn *moldyn,int element1,int element2) { t_albe_mult_params *p; + // set cutoff before parameters (actually only necessary for some pots) + if(moldyn->cutoff==0.0) { + printf("[albe] 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) { diff --git a/potentials/harmonic_oscillator.c b/potentials/harmonic_oscillator.c index b203453..611f518 100644 --- a/potentials/harmonic_oscillator.c +++ b/potentials/harmonic_oscillator.c @@ -19,6 +19,46 @@ #include "../math/math.h" #include "harmonic_oscillator.h" +int harmonic_oscillator_set_params(t_moldyn *moldyn,int element) { + + t_ho_params *p; + + // set cutoff before parameters (actually only necessary for some pots) + if(moldyn->cutoff==0.0) { + printf("[harmonic oscillator] WARNING: no cutoff!\n"); + return -1; + } + + /* alloc mem for potential parameters */ + moldyn->pot_params=malloc(sizeof(t_ho_params)); + if(moldyn->pot_params==NULL) { + perror("[harmonic oscillator] pot params alloc"); + return -1; + } + + /* these are now ho parameters */ + p=moldyn->pot_params; + + switch(element) { + case SI: + /* type: silicon */ + p->spring_constant=HO_SC_SI; + p->equilibrium_distance=HO_ED_SI; + break; + case C: + /* type: carbon */ + p->spring_constant=HO_SC_C; + p->equilibrium_distance=HO_ED_C; + break; + default: + printf("[harmonic oscillator] WARNING: element\n"); + return -1; + } + + + return 0; +} + int harmonic_oscillator(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { t_ho_params *params; @@ -55,3 +95,19 @@ int harmonic_oscillator(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { return 0; } +int harmonic_oscillator_check_2b_bond(t_moldyn *moldyn, + t_atom *ai,t_atom *aj,u8 bc) { + + t_3dvec distance; + double d; + + v3_sub(&distance,&(aj->r),&(ai->r)); + if(bc) check_per_bound(moldyn,&distance); + d=v3_norm(&distance); + + if(d>moldyn->cutoff) + return FALSE; + + return TRUE; +} + diff --git a/potentials/harmonic_oscillator.h b/potentials/harmonic_oscillator.h index c9c9891..13125d9 100644 --- a/potentials/harmonic_oscillator.h +++ b/potentials/harmonic_oscillator.h @@ -15,6 +15,19 @@ typedef struct s_ho_params { } t_ho_params; /* function prototype */ +int harmonic_oscillator_set_params(t_moldyn *moldyn,int element); int harmonic_oscillator(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); +int harmonic_oscillator_check_2b_bond(t_moldyn *moldyn, + t_atom *ai,t_atom *aj,u8 bc); + +/* harmonic oscillator potential parameter defines */ + +// silicon +#define HO_SC_SI 1 +#define HO_ED_SI (0.25*sqrt(3.0)*LC_SI) + +// carbon +#define HO_SC_C 1 +#define HO_ED_C (0.25*sqrt(3.0)*LC_C) #endif diff --git a/potentials/lennard_jones.c b/potentials/lennard_jones.c index 4ff1b43..2424b98 100644 --- a/potentials/lennard_jones.c +++ b/potentials/lennard_jones.c @@ -19,6 +19,62 @@ #include "../math/math.h" #include "lennard_jones.h" +int lennard_jones_set_params(t_moldyn *moldyn,int element) { + + t_lj_params *p; + t_atom a,b; + + // set cutoff before parameters (actually only necessary for some pots) + if(moldyn->cutoff==0.0) { + printf("[lennard jones] WARNING: no cutoff!\n"); + return -1; + } + + /* atoms for correction term */ + a.r.x=0.0; a.r.y=0.0; a.r.z=0.0; + b.r.x=0.0; b.r.y=0.0; b.r.z=moldyn->cutoff; + + /* alloc mem for potential parameters */ + moldyn->pot_params=malloc(sizeof(t_lj_params)); + if(moldyn->pot_params==NULL) { + perror("[lennard jones] pot params alloc"); + return -1; + } + + /* these are now lennard jones parameters */ + p=moldyn->pot_params; + + switch(element) { + case SI: + /* type: silicon */ + p->sigma6=LJ_SIGMA_SI*LJ_SIGMA_SI*LJ_SIGMA_SI; + p->sigma6+=p->sigma6; + p->sigma12=p->sigma6*p->sigma6; + p->epsilon4=4.0*LJ_EPSILON_SI; + p->uc=0.0; // calc it now! + lennard_jones(moldyn,&a,&b,0); + p->uc=moldyn->energy; + moldyn->energy=0.0; + break; + case C: + /* type: carbon */ + p->sigma6=LJ_SIGMA_C*LJ_SIGMA_C*LJ_SIGMA_C; + p->sigma6+=p->sigma6; + p->sigma12=p->sigma6*p->sigma6; + p->epsilon4=4.0*LJ_EPSILON_C; + p->uc=0.0; // calc it now! + lennard_jones(moldyn,&a,&b,0); + p->uc=moldyn->energy; + moldyn->energy=0.0; + break; + default: + printf("[lennard jones] WARNING: element\n"); + return -1; + } + + return 0; +} + int lennard_jones(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { t_lj_params *params; @@ -63,3 +119,17 @@ int lennard_jones(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { return 0; } +int lennard_jones_check_2b_bond(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { + + t_3dvec dist; + double d; + + v3_sub(&dist,&(aj->r),&(ai->r)); + if(bc) check_per_bound(moldyn,&dist); + d=v3_absolute_square(&dist); + + if(d>moldyn->cutoff_square) + return FALSE; + + return TRUE; +} diff --git a/potentials/lennard_jones.h b/potentials/lennard_jones.h index 617658b..8c1c36b 100644 --- a/potentials/lennard_jones.h +++ b/potentials/lennard_jones.h @@ -17,6 +17,18 @@ typedef struct s_lj_params { } t_lj_params; /* function prototype */ +int lennard_jones_set_params(t_moldyn *moldyn,int element); int lennard_jones(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); +int lennard_jones_check_2b_bond(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); + +/* lennard jones potential parameters */ + +// silicon +#define LJ_SIGMA_SI ((0.25*sqrt(3.0)*LC_SI)/1.122462) +#define LJ_EPSILON_SI (2.1678*EV) // TODO + +// carbob +#define LJ_SIGMA_C ((0.25*sqrt(3.0)*LC_C)/1.122462) +#define LJ_EPSILON_C (2.1678*EV) // TODO #endif diff --git a/potentials/tersoff.c b/potentials/tersoff.c index b83bfd6..d322b69 100644 --- a/potentials/tersoff.c +++ b/potentials/tersoff.c @@ -24,6 +24,12 @@ int tersoff_mult_set_params(t_moldyn *moldyn,int element1,int element2) { t_tersoff_mult_params *p; + // set cutoff before parameters (actually only necessary for some pots) + if(moldyn->cutoff==0.0) { + printf("[tersoff] WARNING: no cutoff!\n"); + return -1; + } + /* alloc mem for potential parameters */ moldyn->pot_params=malloc(sizeof(t_tersoff_mult_params)); if(moldyn->pot_params==NULL) { diff --git a/sic.c b/sic.c index 0e14bb1..a3c4f74 100644 --- a/sic.c +++ b/sic.c @@ -262,9 +262,11 @@ int main(int argc,char **argv) { /* choose potential */ #ifdef ALBE - set_potential(&md,MOLDYN_POTENTIAL_AM); + if(set_potential(&md,MOLDYN_POTENTIAL_AM)<0) + return -1; #else - set_potential(&md,MOLDYN_POTENTIAL_TM); + if(set_potential(&md,MOLDYN_POTENTIAL_TM)<0) + return -1; #endif /* cutoff radius & bondlen */ @@ -282,15 +284,17 @@ int main(int argc,char **argv) { * potential parameters */ +#ifndef ALBE /* * tersoff mult potential parameters for SiC */ tersoff_mult_set_params(&md,SI,C); - +#else /* * albe mult potential parameters for SiC */ albe_mult_set_params(&md,SI,C); +#endif /* set (initial) dimensions of simulation volume */ #ifdef ALBE