#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
*/
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)
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->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_k2=albe_mult_3bp_k2;
moldyn->check_2b_bond=albe_mult_check_2b_bond;
break;
+ */
+ case MOLDYN_POTENTIAL_AM:
+ moldyn->func1b=albe_mult_i0;
+ moldyn->func2b=albe_mult_i0_j0;
+ moldyn->func3b_0=albe_mult_i0_j0_k0;
+ moldyn->func3b_1=albe_mult_i0_j1;
+ moldyn->func3b_j1=albe_mult_i0_j2;
+ moldyn->func3b_k2=albe_mult_i0_j2_k0;
+ moldyn->func3b_j2=albe_mult_i0_j3;
+ 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;
t_3dvec orig;
void *ptr;
t_atom *atom;
+ char name[16];
new=a*b*c;
count=moldyn->count;
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);
}
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;ret<new;ret++) {
atom[ret].element=element;
printf("[moldyn] WARNING: forces too high / tau too small!\n");
/* zero absolute time */
- moldyn->time=0.0;
- moldyn->total_steps=0;
+ // should have right values!
+ //moldyn->time=0.0;
+ //moldyn->total_steps=0;
/* debugging, ignore */
moldyn->debug=0;
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);
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,
jtom,
bc_ij);
}
+
+// REWRITE ONCE WORKING!!!
+ /* 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->func3b_0==NULL)
+ continue;
+
+ for(k=0;k<27;k++) {
+
+ bc_ik=(k<dnlc)?0:1;
+#ifdef STATIC_LISTS
+ q=0;
+
+ while(neighbour_i[j][q]!=0) {
+
+ ktom=&(atom[neighbour_i[k][q]]);
+ q++;
+#else
+ that=&(neighbour_i2[k]);
+ list_reset_f(that);
+
+ if(that->start==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->func3b_0(moldyn,
+ &(itom[i]),
+ jtom,
+ ktom,
+ bc_ik|bc_ij);
+#ifdef STATIC_LISTS
+ }
+#else
+ } while(list_next_f(that)!=\
+ L_NO_NEXT_ELEMENT);
+#endif
+
+ }
+
+
+ if(moldyn->func3b_1)
+ moldyn->func3b_1(moldyn,
+ &(itom[i]),
+ jtom,
+ bc_ij);
+
+// UNTIL HERE + change function names!
+
+
+
} while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
#endif
t_list *this;
lc=&(moldyn->lc);
-
- link_cell_init(moldyn,VERBOSE);
-
itom=moldyn->atom;
for(i=0;i<moldyn->count;i++) {
if(moldyn->check_2b_bond(moldyn,itom,jtom,bc)==FALSE)
return 0;
- d=sqrt(d);
-
/* now count this bonding ... */
ba=data;
projects / physik / posic.git / blobdiff