int moldyn_parse_argv(t_moldyn *moldyn,int argc,char **argv) {
int i;
- t_ho_params hop;
- t_lj_params ljp;
- t_tersoff_params tp;
- double s,e;
memset(moldyn,0,sizeof(t_moldyn));
moldyn->tau=MOLDYN_TAU;
moldyn->time_steps=MOLDYN_RUNS;
moldyn->integrate=velocity_verlet;
- moldyn->potential_force_function=lennard_jones;
/* parse argv */
for(i=1;i<argc;i++) {
return 0;
}
-int create_lattice(u8 type,int element,double mass,double lc,
- int a,int b,int c,t_atom **atom) {
+int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
+ u8 attr,u8 bnum,int a,int b,int c) {
int count;
int ret;
t_3dvec origin;
+ t_atom *atom;
count=a*b*c;
+ atom=moldyn->atom;
if(type==FCC) count*=4;
if(type==DIAMOND) count*=8;
- *atom=malloc(count*sizeof(t_atom));
- if(*atom==NULL) {
+ atom=malloc(count*sizeof(t_atom));
+ if(atom==NULL) {
perror("malloc (atoms)");
return -1;
}
switch(type) {
case FCC:
- ret=fcc_init(a,b,c,lc,*atom,&origin);
+ ret=fcc_init(a,b,c,lc,atom,&origin);
break;
case DIAMOND:
- ret=diamond_init(a,b,c,lc,*atom,&origin);
+ ret=diamond_init(a,b,c,lc,atom,&origin);
break;
default:
printf("unknown lattice type (%02x)\n",type);
}
while(count) {
- (*atom)[count-1].element=element;
- (*atom)[count-1].mass=mass;
+ atom[count-1].element=element;
+ atom[count-1].mass=mass;
+ atom[count-1].attr=attr;
+ atom[count-1].bnum=bnum;
count-=1;
}
return ret;
}
-int destroy_lattice(t_atom *atom) {
+int add_atom(t_moldyn *moldyn,int element,double mass,u8 bnum,u8 attr,
+ t_3dvec r,t_3dvec v) {
- if(atom) free(atom);
+ t_atom *atom;
+ void *ptr;
+ int count;
+
+ atom=moldyn->atom;
+ count=++(moldyn->count);
+
+ ptr=realloc(atom,count*sizeof(t_atom));
+ if(!ptr) {
+ perror("[moldyn] realloc (add atom)");
+ return -1;
+ }
+
+ atom=ptr;
+ atom->r=r;
+ atom->v=v;
+ atom->element=element;
+ atom->bnum=bnum;
+ atom->attr=attr;
+
+ return 0;
+}
+
+int destroy_atoms(t_moldyn *moldyn) {
+
+ if(moldyn->atom) free(moldyn->atom);
return 0;
}
return 0;
}
+int moldyn_add_schedule(t_moldyn *moldyn,int runs,double tau ) {
+
+ int count;
+ void *ptr;
+ t_moldyn_schedule *schedule;
+
+ schedule=moldyn->schedule;
+ count=++(schedule->content_count);
+
+ ptr=realloc(moldyn->schedule.runs,count*sizeof(int));
+ if(!ptr) {
+ perror("[moldyn] realloc (runs)");
+ return -1;
+ }
+ moldyn->schedule.runs[count-1]=runs;
+
+ ptr=realloc(schedule->tau,count*sizeof(double));
+ if(!ptr) {
+ perror("[moldyn] realloc (tau)");
+ return -1;
+ }
+ moldyn->schedule.tau[count-1]=tau;
+
+ return 0;
+}
+
+int moldyn_set_schedule_hook(t_moldyn *moldyn,void *hook,void *hook_params) {
+
+ moldyn->schedule.hook=hook;
+ moldyn->schedule.hook_params=hook_params;
+
+ return 0;
+}
+
/*
*
* 'integration of newtons equation' - algorithms
int moldyn_integrate(t_moldyn *moldyn) {
- int i;
+ int i,sched;
unsigned int e,m,s,d,v;
t_3dvec p;
/* calculate initial forces */
moldyn->potential_force_function(moldyn);
+ for(sched=0;sched<moldyn->schedule.content_count;sched++) {
+
+ /* setting amont of runs and finite time step size */
+ moldyn->tau=schedule->tau[sched];
+ moldyn->tau_square=moldyn->tau*moldyn->tau;
+ moldyn->timesteps=schedule->runs[sched];
+
+ /* integration according to schedule */
+
for(i=0;i<moldyn->time_steps;i++) {
/* integration step */
}
}
+ /* check for hooks */
+ if(schedule->hook)
+ schedule->hook(moldyn,schedule->hook_params);
+
return 0;
}
exchange->d=&(params->d[num]);
exchange->h=&(params->h[num]);
+ exchange->betan=pow(*(exchange->beta),*(exchange->n));
exchange->c2=params->c[num]*params->c[num];
exchange->d2=params->d[num]*params->d[num];
exchange->c2d2=exchange->c2/exchange->d2;
t_tersoff_mult_params *params;
t_tersoff_exchange *exchange;
t_3dvec dist_ij,dist_ik,dist_jk;
- t_3dvec db_ij,temp,force;
+ t_3dvec temp,force;
double R,S,s_r;
double d_ij,d_ik,d_jk;
double f_c,df_c,b_ij,f_a,df_a;
- double B,mu;
+ double n,c,d,h,neta,betan,betan_1;
+ double theta,cos_theta,sin_theta;
int num;
params=moldyn->pot_params;
d_ij=exchange->d_ij;
d_ij2=exchange->d_ij2;
- B=*(params->exchange.B);
- mu=*(params->exchange.mu);
-
f_a=params->exchange.f_a;
df_a=params->exchange.df_a;
if(bc) check_per_bound(moldyn,&dist_jk);
d_jk=v3_norm(&dist_jk);
-
- // GO ON HERE !!!
-
- cos_theta=(d_ij2+d_ik*d_ik-d_jk*d_jk)/(2*d_ij*d_ik);
+ beta=*(exchange->beta);
+ betan=exchange->betan;
+ n=*(exchange->n);
+ c=*(exchange->c);
+ d=*(exchange->d);
+ h=*(exchange->h);
+ c2=exchange->c2;
+ d2=exchange->d2;
+ c2d2=exchange->c2d2;
+
+ numer=d_ij2+d_ik*d_ik-d_jk*d_jk;
+ denom=2*d_ij*d_ik;
+ cos_theta=numer/denom;
sin_theta=sqrt(1.0-(cos_theta*cos_theta));
theta=arccos(cos_theta);
+ d_theta=(-1.0/sqrt(1.0-cos_theta*cos_theta))/(denom*denom);
+ d_theta1=2*denom-numer*2*d_ik/d_ij;
+ d_theta2=2*denom-numer*2*d_ij/d_ik;
+ d_theta1*=d_theta;
+ d_theta2*=d_theta;
h_cos=(h-cos_theta);
h_cos2=h_cos*h_cos;
d2_h_cos2=d2-h_cos2;
+ /* some usefull expressions */
+ frac1=c2/(d2-h_cos2);
+ bracket1=1+c2d2-frac1;
+ bracket2=f_c_ik*bracket1;
+ bracket2_n_1=pow(bracket2,n-1.0);
+ bracket2_n=bracket2_n_1*bracket2;
+ bracket3=1+betan*bracket2_n;
+ bracket3_pow_1=pow(bracket3,(-1.0/(2.0*n))-1.0);
+ bracket3_pow=bracket3_pow_1*bracket3;
+
+ /* now go on with calc of b_ij and derivation of b_ij */
+ b_ij=chi*bracket3_pow;
+
+ /* derivation of theta */
+ v3_scale(&force,&dist_ij,d1_theta);
+ v3_scale(&temp,&dist_ik,d_theta2);
+ v3_add(&force,&force,&temp);
+
+ /* part 1 of derivation of b_ij */
+ v3_scale(&force,sin_theta*2*h_cos*f_c_ik*frac1);
+
+ /* part 2 of derivation of b_ij */
+ v3_scale(&temp,&dist_ik,df_c_ik*bracket1);
+
+ /* sum up and scale ... */
+ v3_add(&temp,&temp,&force);
+ scale=bracket2_n_1*n*betan*(1+betan*bracket3_pow_1)*chi*(1.0/(2.0*n));
+ v3_scale(&temp,&temp,scale);
+
+ /* now construct an energy and a force out of that */
+ v3_scale(&temp,&temp,f_a);
+ v3_scale(&force,&dist_ij,df_a*b_ij);
+ v3_add(&temp,&temp,&force);
+ v3_scale(&temp,&temp,f_c);
+ v3_scale(&force,&dist_ij,df_c*b_ij*f_a);
+ v3_add(&force,&force,&temp);
+
/* add forces */
v3_add(&(ai->f),&(ai->f),&force);
/* energy is 0.5 f_r f_c, but we will sum it up twice ... */
projects / physik / posic.git / blobdiff