X-Git-Url: https://hackdaworld.org/gitweb/?a=blobdiff_plain;f=moldyn.c;h=1aff5a8df207c2b5d22077c70ccbc31c128afdd6;hb=2b6a06cdf69b0d8f1dafc0e21b89c6bb4bfc192c;hp=988b4671780d6fb68a8826c92390c2cedd2ff663;hpb=d7f67c88195ab155f2737e57cc5e81973d3feb0c;p=physik%2Fposic.git diff --git a/moldyn.c b/moldyn.c index 988b467..1aff5a8 100644 --- a/moldyn.c +++ b/moldyn.c @@ -30,99 +30,6 @@ #include "potentials/tersoff.h" #endif - -/* - * global variables, pse and atom colors (only needed here) - */ - -static char *pse_name[]={ - "*", - "H", - "He", - "Li", - "Be", - "B", - "C", - "N", - "O", - "F", - "Ne", - "Na", - "Mg", - "Al", - "Si", - "P", - "S", - "Cl", - "Ar", -}; - -static char *pse_col[]={ - "*", - "White", - "He", - "Li", - "Be", - "B", - "Gray", - "N", - "Blue", - "F", - "Ne", - "Na", - "Mg", - "Al", - "Yellow", - "P", - "S", - "Cl", - "Ar", -}; - -static double *pse_mass[]={ - 0, - 0, - 0, - 0, - 0, - 0, - M_C, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - M_SI, - 0, - 0, - 0, - 0, -}; - -static double *pse_lc[]={ - 0, - 0, - 0, - 0, - 0, - 0, - LC_C, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - LC_SI, - 0, - 0, - 0, - 0, -}; - /* * the moldyn functions */ @@ -202,37 +109,27 @@ int set_pressure(t_moldyn *moldyn,double p_ref) { int set_p_scale(t_moldyn *moldyn,u8 ptype,double ptc) { - moldyn->pt_scalei&=(^(P_SCALE_MASK)); + moldyn->pt_scale&=(~(P_SCALE_MASK)); moldyn->pt_scale|=ptype; moldyn->p_tc=ptc; - printf("[moldyn] p/t scaling:\n"); + printf("[moldyn] p scaling:\n"); printf(" p: %s",ptype?"yes":"no "); if(ptype) printf(" | type: %02x | factor: %f",ptype,ptc); printf("\n"); - printf(" t: %s",ttype?"yes":"no "); - if(ttype) - printf(" | type: %02x | factor: %f",ttype,ttc); - printf("\n"); - return 0; } int set_t_scale(t_moldyn *moldyn,u8 ttype,double ttc) { - moldyn->pt_scalei&=(^(T_SCALE_MASK)); + moldyn->pt_scale&=(~(T_SCALE_MASK)); moldyn->pt_scale|=ttype; moldyn->t_tc=ttc; - printf("[moldyn] p/t scaling:\n"); - - printf(" p: %s",ptype?"yes":"no "); - if(ptype) - printf(" | type: %02x | factor: %f",ptype,ptc); - printf("\n"); + printf("[moldyn] t scaling:\n"); printf(" t: %s",ttype?"yes":"no "); if(ttype) @@ -318,13 +215,14 @@ int set_potential(t_moldyn *moldyn,u8 type) { switch(type) { case MOLDYN_POTENTIAL_TM: - moldyn->func1b=tersoff_mult_1bp; - moldyn->func3b_j1=tersoff_mult_3bp_j1; - moldyn->func3b_k1=tersoff_mult_3bp_k1; - moldyn->func3b_j2=tersoff_mult_3bp_j2; - moldyn->func3b_k2=tersoff_mult_3bp_k2; - // missing: check 2b bond func + moldyn->func_i0=tersoff_mult_1bp; + moldyn->func_j1=tersoff_mult_3bp_j1; + moldyn->func_j1_k0=tersoff_mult_3bp_k1; + moldyn->func_j1c=tersoff_mult_3bp_j2; + 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; @@ -332,12 +230,23 @@ 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_AM: + moldyn->func_i0=albe_mult_i0; + moldyn->func_j0=albe_mult_i0_j0; + moldyn->func_j0_k0=albe_mult_i0_j0_k0; + moldyn->func_j0e=albe_mult_i0_j1; + moldyn->func_j1=albe_mult_i0_j2; + moldyn->func_j1_k0=albe_mult_i0_j2_k0; + moldyn->func_j1c=albe_mult_i0_j3; + moldyn->check_2b_bond=albe_mult_check_2b_bond; + break; case MOLDYN_POTENTIAL_HO: - moldyn->func2b=harmonic_oscillator; + moldyn->func_j0=harmonic_oscillator; moldyn->check_2b_bond=harmonic_oscillator_check_2b_bond; break; case MOLDYN_POTENTIAL_LJ: - moldyn->func2b=lennard_jones; + moldyn->func_j0=lennard_jones; moldyn->check_2b_bond=lennard_jones_check_2b_bond; break; default: @@ -582,6 +491,7 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass, t_3dvec orig; void *ptr; t_atom *atom; + char name[16]; new=a*b*c; count=moldyn->count; @@ -616,18 +526,21 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass, case CUBIC: set_nn_dist(moldyn,lc); ret=cubic_init(a,b,c,lc,atom,&orig); + strcpy(name,"cubic"); break; case FCC: 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,&orig); + strcpy(name,"fcc"); break; case DIAMOND: 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,&orig); + strcpy(name,"diamond"); break; default: printf("unknown lattice type (%02x)\n",type); @@ -644,7 +557,7 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass, } moldyn->count+=new; - printf("[moldyn] created lattice with %d atoms\n",new); + printf("[moldyn] created %s lattice with %d atoms\n",name,new); for(ret=0;rettime=0.0; - moldyn->total_steps=0; + // should have right values! + //moldyn->time=0.0; + //moldyn->total_steps=0; /* debugging, ignore */ moldyn->debug=0; @@ -1711,7 +1625,11 @@ int moldyn_integrate(t_moldyn *moldyn) { temperature_calc(moldyn); virial_sum(moldyn); pressure_calc(moldyn); - //thermodynamic_pressure_calc(moldyn); + /* + thermodynamic_pressure_calc(moldyn); + printf("\n\nDEBUG: numeric pressure calc: %f\n\n", + moldyn->tp/BAR); + */ /* calculate fluctuations + averages */ average_and_fluctuation_calc(moldyn); @@ -1942,8 +1860,8 @@ int potential_force_calc(t_moldyn *moldyn) { /* single particle potential/force */ if(itom[i].attr&ATOM_ATTR_1BP) - if(moldyn->func1b) - moldyn->func1b(moldyn,&(itom[i])); + if(moldyn->func_i0) + moldyn->func_i0(moldyn,&(itom[i])); if(!(itom[i].attr&(ATOM_ATTR_2BP|ATOM_ATTR_3BP))) continue; @@ -1958,8 +1876,14 @@ int potential_force_calc(t_moldyn *moldyn) { dnlc=lc->dnlc; +#ifndef STATIC_LISTS + /* check for later 3 body interaction */ + if(itom[i].attr&ATOM_ATTR_3BP) + memcpy(neighbour_i2,neighbour_i,27*sizeof(t_list)); +#endif + /* first loop over atoms j */ - if(moldyn->func2b) { + if(moldyn->func_j0) { for(j=0;j<27;j++) { bc_ij=(jattr&ATOM_ATTR_2BP)& - (itom[i].attr&ATOM_ATTR_2BP)) { - moldyn->func2b(moldyn, - &(itom[i]), - jtom, - bc_ij); - } - } #else this=&(neighbour_i[j]); list_reset_f(this); @@ -1991,35 +1903,103 @@ int potential_force_calc(t_moldyn *moldyn) { do { jtom=this->current->data; +#endif if(jtom==&(itom[i])) continue; + /* reset 3bp run */ + moldyn->run3bp=1; + if((jtom->attr&ATOM_ATTR_2BP)& (itom[i].attr&ATOM_ATTR_2BP)) { - moldyn->func2b(moldyn, - &(itom[i]), - jtom, - bc_ij); + moldyn->func_j0(moldyn, + &(itom[i]), + jtom, + bc_ij); + } + + /* 3 body potential/force */ + + /* in j loop, 3bp run can be skipped */ + if(!(moldyn->run3bp)) + continue; + + if(!(itom[i].attr&ATOM_ATTR_3BP)) + continue; + + if(!(jtom->attr&ATOM_ATTR_3BP)) + continue; + + if(moldyn->func_j0_k0==NULL) + continue; + + /* first loop over atoms k in first j loop */ + for(k=0;k<27;k++) { + + bc_ik=(kstart==NULL) + continue; + + do { + ktom=that->current->data; +#endif + + if(!(ktom->attr&ATOM_ATTR_3BP)) + continue; + + //if(ktom==jtom) + // continue; + + if(ktom==&(itom[i])) + continue; + + moldyn->func_j0_k0(moldyn, + &(itom[i]), + jtom, + ktom, + bc_ik|bc_ij); +#ifdef STATIC_LISTS } +#else + } while(list_next_f(that)!=\ + L_NO_NEXT_ELEMENT); +#endif + + } + + /* finish of first j loop after first k loop */ + if(moldyn->func_j0e) + moldyn->func_j0e(moldyn, + &(itom[i]), + jtom, + bc_ij); + +#ifdef STATIC_LISTS + } +#else } while(list_next_f(this)!=L_NO_NEXT_ELEMENT); #endif } } - /* 3 body potential/force */ + /* continued 3 body potential/force */ if(!(itom[i].attr&ATOM_ATTR_3BP)) continue; - /* copy the neighbour lists */ -#ifdef STATIC_LISTS - /* no copy needed for static lists */ -#else - memcpy(neighbour_i2,neighbour_i,27*sizeof(t_list)); -#endif - /* second loop over atoms j */ for(j=0;j<27;j++) { @@ -2052,18 +2032,18 @@ int potential_force_calc(t_moldyn *moldyn) { /* reset 3bp run */ moldyn->run3bp=1; - if(moldyn->func3b_j1) - moldyn->func3b_j1(moldyn, - &(itom[i]), - jtom, - bc_ij); + if(moldyn->func_j1) + moldyn->func_j1(moldyn, + &(itom[i]), + jtom, + bc_ij); - /* in first j loop, 3bp run can be skipped */ + /* in j loop, 3bp run can be skipped */ if(!(moldyn->run3bp)) continue; - /* first loop over atoms k */ - if(moldyn->func3b_k1) { + /* first loop over atoms k in second j loop */ + if(moldyn->func_j1_k0) { for(k=0;k<27;k++) { @@ -2089,17 +2069,17 @@ int potential_force_calc(t_moldyn *moldyn) { if(!(ktom->attr&ATOM_ATTR_3BP)) continue; - if(ktom==jtom) - continue; + //if(ktom==jtom) + // continue; if(ktom==&(itom[i])) continue; - moldyn->func3b_k1(moldyn, - &(itom[i]), - jtom, - ktom, - bc_ik|bc_ij); + moldyn->func_j1_k0(moldyn, + &(itom[i]), + jtom, + ktom, + bc_ik|bc_ij); #ifdef STATIC_LISTS } #else @@ -2111,14 +2091,15 @@ int potential_force_calc(t_moldyn *moldyn) { } - if(moldyn->func3b_j2) - moldyn->func3b_j2(moldyn, - &(itom[i]), - jtom, - bc_ij); + /* continued j loop after first k loop */ + if(moldyn->func_j1c) + moldyn->func_j1c(moldyn, + &(itom[i]), + jtom, + bc_ij); /* second loop over atoms k */ - if(moldyn->func3b_k2) { + if(moldyn->func_j1_k1) { for(k=0;k<27;k++) { @@ -2144,17 +2125,17 @@ int potential_force_calc(t_moldyn *moldyn) { if(!(ktom->attr&ATOM_ATTR_3BP)) continue; - if(ktom==jtom) - continue; + //if(ktom==jtom) + // continue; if(ktom==&(itom[i])) continue; - moldyn->func3b_k2(moldyn, - &(itom[i]), - jtom, - ktom, - bc_ik|bc_ij); + moldyn->func_j1_k1(moldyn, + &(itom[i]), + jtom, + ktom, + bc_ik|bc_ij); #ifdef STATIC_LISTS } @@ -2167,11 +2148,11 @@ int potential_force_calc(t_moldyn *moldyn) { } - /* 2bp post function */ - if(moldyn->func3b_j3) { - moldyn->func3b_j3(moldyn, - &(itom[i]), - jtom,bc_ij); + /* finish of j loop after second k loop */ + if(moldyn->func_j1e) { + moldyn->func_j1e(moldyn, + &(itom[i]), + jtom,bc_ij); } #ifdef STATIC_LISTS } @@ -2180,7 +2161,7 @@ int potential_force_calc(t_moldyn *moldyn) { #endif } - + #ifdef DEBUG //printf("\n\n"); #endif @@ -2422,9 +2403,6 @@ int process_2b_bonds(t_moldyn *moldyn,void *data, t_list *this; lc=&(moldyn->lc); - - link_cell_init(moldyn,VERBOSE); - itom=moldyn->atom; for(i=0;icount;i++) { @@ -2684,8 +2662,6 @@ int bond_analyze_process(t_moldyn *moldyn,t_atom *itom,t_atom *jtom, if(moldyn->check_2b_bond(moldyn,itom,jtom,bc)==FALSE) return 0; - d=sqrt(d); - /* now count this bonding ... */ ba=data;