X-Git-Url: https://hackdaworld.org/gitweb/?p=physik%2Fposic.git;a=blobdiff_plain;f=moldyn.c;h=1d6e0b7dec9e3032df43f33baa692ba078bb5f7a;hp=d715c52a11f91aa9f43244b014c281f3c2a00cc4;hb=HEAD;hpb=fe84cdc3fc897581636aff89b6f4fd12db1d2718 diff --git a/moldyn.c b/moldyn.c index d715c52..1d6e0b7 100644 --- a/moldyn.c +++ b/moldyn.c @@ -34,6 +34,7 @@ #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 @@ -270,19 +271,29 @@ int set_potential(t_moldyn *moldyn,u8 type) { 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: @@ -2319,6 +2330,7 @@ int velocity_verlet(t_moldyn *moldyn) { /* constraint relaxation */ if(crtt) { + // forces basis_trafo(&(atom[i].f),FORWARD, trafo_angle[2*i],trafo_angle[2*i+1]); if(constraints[3*i]) @@ -2329,6 +2341,17 @@ int velocity_verlet(t_moldyn *moldyn) { 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 @@ -2369,6 +2392,7 @@ int velocity_verlet(t_moldyn *moldyn) { /* constraint relaxation */ if(crtt) { + // forces basis_trafo(&(atom[i].f),FORWARD, trafo_angle[2*i],trafo_angle[2*i+1]); if(constraints[3*i]) @@ -2379,6 +2403,17 @@ int velocity_verlet(t_moldyn *moldyn) { 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)] */ @@ -2463,8 +2498,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; @@ -2479,8 +2514,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=(jrun3bp=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 } #ifdef STATIC_LISTS } @@ -2530,7 +2628,7 @@ int potential_force_calc(t_moldyn *moldyn) { } } - /* 3 body potential/force */ + /* continued 3 body potential/force */ if(!(itom[i].attr&ATOM_ATTR_3BP)) continue; @@ -2583,18 +2681,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++) { @@ -2627,8 +2725,8 @@ 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; @@ -2652,14 +2750,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++) { @@ -2692,17 +2791,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 } @@ -2717,11 +2816,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 } @@ -2732,7 +2831,7 @@ int potential_force_calc(t_moldyn *moldyn) { #endif } - + #ifdef DEBUG //printf("\n\n"); #endif