#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
pthread_mutex_t emutex;
#endif
+/* fully constrained relaxation technique - global pointers */
+u8 crtt;
+u8 *constraints;
+double *trafo_angle;
+
/*
* the moldyn functions
*/
switch(type) {
case MOLDYN_POTENTIAL_TM:
- moldyn->func1b=tersoff_mult_1bp;
+ //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;
moldyn->check_2b_bond=tersoff_mult_check_2b_bond;
break;
+ 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->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->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:
int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,
u8 attr,u8 brand,int a,int b,int c,t_3dvec *origin,
- t_part_params *p_params,t_defect_params *d_params) {
+ t_part_params *p_params,t_defect_params *d_params,
+ t_offset_params *o_params) {
int new,count;
int ret;
switch(type) {
case CUBIC:
+ if(o_params->use)
+ v3_add(&orig,&orig,&(o_params->o));
set_nn_dist(moldyn,lc);
ret=cubic_init(a,b,c,lc,atom,&orig,p_params,d_params);
strcpy(name,"cubic");
case FCC:
if(!origin)
v3_scale(&orig,&orig,0.5);
+ if(o_params->use)
+ v3_add(&orig,&orig,&(o_params->o));
set_nn_dist(moldyn,0.5*sqrt(2.0)*lc);
ret=fcc_init(a,b,c,lc,atom,&orig,p_params,d_params);
strcpy(name,"fcc");
case DIAMOND:
if(!origin)
v3_scale(&orig,&orig,0.25);
+ if(o_params->use)
+ v3_add(&orig,&orig,&(o_params->o));
set_nn_dist(moldyn,0.25*sqrt(3.0)*lc);
ret=diamond_init(a,b,c,lc,atom,&orig,p_params,d_params);
strcpy(name,"diamond");
t_atom *new,*old;
int cnt;
+#if defined LOWMEM_LISTS || defined PTHREADS
+ void *ptr;
+#endif
old=moldyn->atom;
free(old);
+#ifdef LOWMEM_LISTS
+ ptr=realloc(moldyn->lc.subcell->list,moldyn->count*sizeof(int));
+ if(!ptr) {
+ perror("[moldyn] list realloc (del atom)");
+ return -1;
+ }
+ moldyn->lc.subcell->list=ptr;
+ // update lists
+ link_cell_update(moldyn);
+#endif
+
+#ifdef PTHREADS
+ ptr=realloc(amutex,moldyn->count*sizeof(pthread_mutex_t));
+ if(!ptr) {
+ perror("[moldyn] mutex realloc (add atom)");
+ return -1;
+ }
+ amutex=ptr;
+ pthread_mutex_destroy(&(amutex[moldyn->count+1]));
+#endif
+
+
return 0;
}
temperature_calc(moldyn);
virial_sum(moldyn);
pressure_calc(moldyn);
- /*
+#ifdef PDEBUG
thermodynamic_pressure_calc(moldyn);
printf("\n\nDEBUG: numeric pressure calc: %f\n\n",
moldyn->tp/BAR);
- */
+#endif
/* calculate fluctuations + averages */
average_and_fluctuation_calc(moldyn);
if(e) {
if(!(moldyn->total_steps%e))
dprintf(moldyn->efd,
- "%f %f %f %f\n",
+ "%f %f %f %f %f %f\n",
moldyn->time,moldyn->ekin/energy_scale,
moldyn->energy/energy_scale,
- get_total_energy(moldyn)/energy_scale);
+ get_total_energy(moldyn)/energy_scale,
+ moldyn->ekin/EV,moldyn->energy/EV);
}
if(m) {
if(!(moldyn->total_steps%m)) {
}
/* display progress */
+#ifndef PDEBUG
if(!(i%10)) {
+#endif
/* get current time */
gettimeofday(&t2,NULL);
printf("sched:%d, steps:%d/%d, T:%4.1f/%4.1f P:%4.1f/%4.1f V:%6.1f (%d)\n",
sched->count,i,moldyn->total_steps,
moldyn->t,moldyn->t_avg,
+#ifndef PDEBUG
moldyn->p/BAR,moldyn->p_avg/BAR,
- //moldyn->p/BAR,(moldyn->p-2.0*moldyn->ekin/(3.0*moldyn->volume))/BAR,
+#else
+ moldyn->p/BAR,(moldyn->p-2.0*moldyn->ekin/(3.0*moldyn->volume))/BAR,
+#endif
moldyn->volume,
(int)(t2.tv_sec-t1.tv_sec));
/* copy over time */
t1=t2;
+#ifndef PDEBUG
}
+#endif
/* increase absolute time */
moldyn->time+=moldyn->tau;
}
+ /* writing a final save file! */
+ if(s) {
+ snprintf(dir,128,"%s/s-final.save",moldyn->vlsdir);
+ fd=open(dir,O_WRONLY|O_TRUNC|O_CREAT,S_IRUSR|S_IWUSR);
+ if(fd<0) perror("[moldyn] save fd open");
+ else {
+ write(fd,moldyn,sizeof(t_moldyn));
+ write(fd,moldyn->atom,
+ moldyn->count*sizeof(t_atom));
+ }
+ close(fd);
+ }
+ /* writing a final visual file! */
+ if(a)
+ visual_atoms(moldyn);
+
+ return 0;
+}
+
+/* basis trafo */
+
+#define FORWARD 0
+#define BACKWARD 1
+
+int basis_trafo(t_3dvec *r,u8 dir,double z,double x) {
+
+ t_3dvec tmp;
+
+ if(dir==FORWARD) {
+ if(z!=0.0) {
+ v3_copy(&tmp,r);
+ r->x=cos(z)*tmp.x-sin(z)*tmp.y;
+ r->y=sin(z)*tmp.x+cos(z)*tmp.y;
+ }
+ if(x!=0.0) {
+ v3_copy(&tmp,r);
+ r->y=cos(x)*tmp.y-sin(x)*tmp.z;
+ r->z=sin(x)*tmp.y+cos(x)*tmp.z;
+ }
+ }
+ else {
+ if(x!=0.0) {
+ v3_copy(&tmp,r);
+ r->y=cos(-x)*tmp.y-sin(-x)*tmp.z;
+ r->z=sin(-x)*tmp.y+cos(-x)*tmp.z;
+ }
+ if(z!=0.0) {
+ v3_copy(&tmp,r);
+ r->x=cos(-z)*tmp.x-sin(-z)*tmp.y;
+ r->y=sin(-z)*tmp.x+cos(-z)*tmp.y;
+ }
+ }
+
return 0;
}
tau_square=moldyn->tau_square;
for(i=0;i<count;i++) {
+
/* check whether fixed atom */
if(atom[i].attr&ATOM_ATTR_FP)
continue;
+
/* new positions */
h=0.5/atom[i].mass;
+
+ /* constraint relaxation */
+ if(crtt) {
+ // forces
+ basis_trafo(&(atom[i].f),FORWARD,
+ trafo_angle[2*i],trafo_angle[2*i+1]);
+ if(constraints[3*i])
+ atom[i].f.x=0;
+ if(constraints[3*i+1])
+ atom[i].f.y=0;
+ if(constraints[3*i+2])
+ atom[i].f.z=0;
+ basis_trafo(&(atom[i].f),BACKWARD,
+ trafo_angle[2*i],trafo_angle[2*i+1]);
+ // velocities
+ basis_trafo(&(atom[i].v),FORWARD,
+ trafo_angle[2*i],trafo_angle[2*i+1]);
+ if(constraints[3*i])
+ atom[i].v.x=0;
+ if(constraints[3*i+1])
+ atom[i].v.y=0;
+ if(constraints[3*i+2])
+ atom[i].v.z=0;
+ basis_trafo(&(atom[i].v),BACKWARD,
+ trafo_angle[2*i],trafo_angle[2*i+1]);
+ }
+
+#ifndef QUENCH
v3_scale(&delta,&(atom[i].v),tau);
v3_add(&(atom[i].r),&(atom[i].r),&delta);
+#endif
v3_scale(&delta,&(atom[i].f),h*tau_square);
v3_add(&(atom[i].r),&(atom[i].r),&delta);
+ //check_per_bound_and_care_for_pbc(moldyn,&(atom[i]));
check_per_bound(moldyn,&(atom[i].r));
/* velocities [actually v(t+tau/2)] */
/* forces depending on chosen potential */
#ifndef ALBE_FAST
- potential_force_calc(moldyn);
+ // if albe, use albe force calc routine
+ //if(moldyn->func3b_j1==albe_mult_3bp_j1)
+ // albe_potential_force_calc(moldyn);
+ //else
+ potential_force_calc(moldyn);
#else
albe_potential_force_calc(moldyn);
#endif
/* check whether fixed atom */
if(atom[i].attr&ATOM_ATTR_FP)
continue;
+
+ /* constraint relaxation */
+ if(crtt) {
+ // forces
+ basis_trafo(&(atom[i].f),FORWARD,
+ trafo_angle[2*i],trafo_angle[2*i+1]);
+ if(constraints[3*i])
+ atom[i].f.x=0;
+ if(constraints[3*i+1])
+ atom[i].f.y=0;
+ if(constraints[3*i+2])
+ atom[i].f.z=0;
+ basis_trafo(&(atom[i].f),BACKWARD,
+ trafo_angle[2*i],trafo_angle[2*i+1]);
+ // velocities
+ basis_trafo(&(atom[i].v),FORWARD,
+ trafo_angle[2*i],trafo_angle[2*i+1]);
+ if(constraints[3*i])
+ atom[i].v.x=0;
+ if(constraints[3*i+1])
+ atom[i].v.y=0;
+ if(constraints[3*i+2])
+ atom[i].v.z=0;
+ basis_trafo(&(atom[i].v),BACKWARD,
+ trafo_angle[2*i],trafo_angle[2*i+1]);
+ }
+
/* again velocities [actually v(t+tau)] */
v3_scale(&delta,&(atom[i].f),0.5*tau/atom[i].mass);
v3_add(&(atom[i].v),&(atom[i].v),&delta);
/* 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;
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=(j<dnlc)?0:1;
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=(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->func_j0_k0(moldyn,
+ &(itom[i]),
+ jtom,
+ ktom,
+ bc_ik|bc_ij);
+#ifdef STATIC_LISTS
}
#ifdef STATIC_LISTS
}
}
}
- /* 3 body potential/force */
+ /* continued 3 body potential/force */
if(!(itom[i].attr&ATOM_ATTR_3BP))
continue;
/* 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++) {
if(!(ktom->attr&ATOM_ATTR_3BP))
continue;
- if(ktom==jtom)
- continue;
+ //if(ktom==jtom)
+ // continue;
if(ktom==&(itom[i]))
continue;
}
- 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++) {
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
}
}
- /* 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
}
#endif
}
-
+
#ifdef DEBUG
//printf("\n\n");
#endif
return 0;
}
+int check_per_bound_and_care_for_pbc(t_moldyn *moldyn,t_atom *a) {
+
+ double x,y,z;
+ t_3dvec *dim;
+
+ dim=&(moldyn->dim);
+
+ x=dim->x/2;
+ y=dim->y/2;
+ z=dim->z/2;
+
+ if(moldyn->status&MOLDYN_STAT_PBX) {
+ if(a->r.x>=x) {
+ a->pbc[0]+=1;
+ a->r.x-=dim->x;
+ }
+ else if(-a->r.x>x) {
+ a->pbc[0]-=1;
+ a->r.x+=dim->x;
+ }
+ }
+ if(moldyn->status&MOLDYN_STAT_PBY) {
+ if(a->r.y>=y) {
+ a->pbc[1]+=1;
+ a->r.y-=dim->y;
+ }
+ else if(-a->r.y>y) {
+ a->pbc[1]-=1;
+ a->r.y+=dim->y;
+ }
+ }
+ if(moldyn->status&MOLDYN_STAT_PBZ) {
+ if(a->r.z>=z) {
+ a->pbc[2]+=1;
+ a->r.z-=dim->z;
+ }
+ else if(-a->r.z>z) {
+ a->pbc[2]-=1;
+ a->r.z+=dim->z;
+ }
+ }
+
+ return 0;
+}
+
/*
* debugging / critical check functions
*/
if(fsize!=sizeof(t_moldyn)+size) {
corr=fsize-sizeof(t_moldyn)-size;
printf("[moldyn] WARNING: lsf (illegal file size)\n");
- printf(" moifying offset:\n");
+ printf(" modifying offset:\n");
printf(" - current pos: %d\n",sizeof(t_moldyn));
printf(" - atom size: %d\n",size);
printf(" - file size: %d\n",fsize);
int i;
t_atom *atom;
t_3dvec dist;
+ t_3dvec final_r;
double d2;
int a_cnt;
int b_cnt;
for(i=0;i<moldyn->count;i++) {
- v3_sub(&dist,&(atom[i].r),&(atom[i].r_0));
- check_per_bound(moldyn,&dist);
+ /* care for pb crossing */
+ final_r.x=atom[i].r.x+atom[i].pbc[0]*moldyn->dim.x;
+ final_r.y=atom[i].r.y+atom[i].pbc[1]*moldyn->dim.y;
+ final_r.z=atom[i].r.z+atom[i].pbc[2]*moldyn->dim.z;
+ /* calculate distance */
+ v3_sub(&dist,&final_r,&(atom[i].r_0));
d2=v3_absolute_square(&dist);
if(atom[i].brand) {
return 0;
}
+int calculate_msd(t_moldyn *moldyn,double *msd) {
+
+ int i;
+ t_atom *atom;
+ t_3dvec dist;
+ t_3dvec final_r;
+ double d2;
+ int a_cnt;
+ int b_cnt;
+
+ atom=moldyn->atom;
+ msd[0]=0;
+ msd[1]=0;
+ msd[2]=0;
+ a_cnt=0;
+ b_cnt=0;
+
+ for(i=0;i<moldyn->count;i++) {
+
+ /* care for pb crossing */
+ if(atom[i].pbc[0]|atom[i].pbc[1]|atom[i].pbc[2]) {
+ printf("[moldyn] msd pb crossings for atom %d\n",i);
+ printf(" x: %d y: %d z: %d\n",
+ atom[i].pbc[0],atom[i].pbc[1],atom[i].pbc[2]);
+ }
+ final_r.x=atom[i].r.x+atom[i].pbc[0]*moldyn->dim.x;
+ final_r.y=atom[i].r.y+atom[i].pbc[1]*moldyn->dim.y;
+ final_r.z=atom[i].r.z+atom[i].pbc[2]*moldyn->dim.z;
+ /* calculate distance */
+ v3_sub(&dist,&final_r,&(atom[i].r_0));
+ d2=v3_absolute_square(&dist);
+
+ if(atom[i].brand) {
+ b_cnt+=1;
+ msd[1]+=d2;
+ }
+ else {
+ a_cnt+=1;
+ msd[0]+=d2;
+ }
+
+ msd[2]+=d2;
+ }
+
+ msd[0]/=a_cnt;
+ msd[1]/=b_cnt;
+ msd[2]/=moldyn->count;
+
+ return 0;
+}
+
int bonding_analyze(t_moldyn *moldyn,double *cnt) {
return 0;
ba=data;
/* increase total bond counter
- * ... double counting!
*/
- ba->tcnt+=2;
+ ba->tcnt+=1;
if(itom->brand==0)
ba->acnt[jtom->tag]+=1;
int bond_analyze(t_moldyn *moldyn,double *quality) {
- // by now: # bonds of type 'a-4b' and 'b-4a' / # bonds total
-
- int qcnt;
- int ccnt,cset;
+ int qcnt4;
+ int qcnt3;
+ int ccnt4;
+ int ccnt3;
+ int bcnt;
t_ba ba;
int i;
t_atom *atom;
memset(ba.bcnt,0,moldyn->count*sizeof(int));
ba.tcnt=0;
- qcnt=0;
- ccnt=0;
- cset=0;
+ qcnt4=0; qcnt3=0;
+ ccnt4=0; ccnt3=0;
+ bcnt=0;
atom=moldyn->atom;
for(i=0;i<moldyn->count;i++) {
if(atom[i].brand==0) {
if((ba.acnt[i]==0)&(ba.bcnt[i]==4))
- qcnt+=4;
+ qcnt4+=4;
+ if((ba.acnt[i]==0)&(ba.bcnt[i]==3))
+ qcnt3+=3;
}
else {
if((ba.acnt[i]==4)&(ba.bcnt[i]==0)) {
- qcnt+=4;
- ccnt+=1;
+ qcnt4+=4;
+ ccnt4+=1;
}
- cset+=1;
+ if((ba.acnt[i]==3)&(ba.bcnt[i]==0)) {
+ qcnt3+=4;
+ ccnt3+=1;
+ }
+ bcnt+=1;
}
}
- printf("[moldyn] bond analyze: c_cnt=%d | set=%d\n",ccnt,cset);
- printf("[moldyn] bond analyze: q_cnt=%d | tot=%d\n",qcnt,ba.tcnt);
-
if(quality) {
- quality[0]=1.0*ccnt/cset;
- quality[1]=1.0*qcnt/ba.tcnt;
+ quality[0]=1.0*ccnt4/bcnt;
+ quality[1]=1.0*ccnt3/bcnt;
}
else {
- printf("[moldyn] bond analyze: c_bnd_q=%f\n",1.0*qcnt/ba.tcnt);
- printf("[moldyn] bond analyze: tot_q=%f\n",1.0*qcnt/ba.tcnt);
+ printf("[moldyn] bond analyze: %f %f\n",
+ 1.0*ccnt4/bcnt,1.0*ccnt3/bcnt);
}
return 0;
projects / physik / posic.git / blobdiff