diff --git a/doc/ug/setup.tex b/doc/ug/setup.tex index c90bbefd9ce..d46266da6ce 100755 --- a/doc/ug/setup.tex +++ b/doc/ug/setup.tex @@ -679,8 +679,10 @@ \section{Creating bonds when particles collide} \begin{essyntax} \variant{1} on\_collision \variant{2} on\_collision off -\variant{4} on\_collision \opt{exception} bind_centers \var{d} \var{bond1} -\variant{5} on\_collision \opt{exception} bind_at_point_of_collision \var{d} \var{bond1} \var{bond2} \var{type} +\variant{3} on\_collision \opt{exception} bind_centers \var{d} \var{bond1} +\variant{4} on\_collision \opt{exception} bind_at_point_of_collision \var{d} \var{bond1} \var{bond2} \var{type} +\variant{5} on\_collision \opt{exception} glue_to_surface \var{d} \var{bond1} \var{bond2} \var{type} \var{type2} \var{type3} \var{type4} \var{d2} +\variant{6} on\_collision \opt{exception} bind_three_particles \var{d} \var{bond1} \var{bond2} \var{res} \end{essyntax} With the help of the feature \feature{COLLISION_DETECTION}, bonds @@ -728,6 +730,17 @@ \section{Creating bonds when particles collide} created at the point of collision (if applicable). Be sure not to define a short-ranged interaction for this particle type, as two particles will be generated in the same place. + +\item \lit{glue_to_surface} is used to fix a particle of type \var{type2} onto the surface of a particle of type \var{type3}. +This is achieved by creating a virtual site (particle type \var{type}) which +is rigidly connected to the particle of \var{type3}. A bond of type +\var{bond2} is then created between the virtual site and the particle +of \var{type}. +Additionally, a bond of type \var{bond1} between the colliding +particles is also created. +After the collision, the particle of type \var{type3} is changed to type \var{type4}. +\item \lit{bind_three_particles} allows for the creating of agglomerates which maintain their shape similarly to those create by the \lit{bind_ad_point_of_collision} method. The present approach works without virtual sites. Instead, for each two-particle collision, the surrounding is searched for a third particle. If one is found, angular bounds are placed on each of the three particles in addition to the distance based bounds between the particle centers. The id of the angular bonds is determined from the angle between the particles. Zero degrees corresponds to bond id \var{bond2}, where as 180 degrees corresponds to bond id \var{bond2} +\var{res}. +this method das not depend on the particles' rotational degrees of freedom being integrated. Virtual sites are also not required, and the method is implemented to run on more than one cpu core. \end{itemize} The code can throw an exception (background error) in case two diff --git a/src/core/actor/DipolarDirectSum.cpp b/src/core/actor/DipolarDirectSum.cpp index 186df140744..a446188b34b 100644 --- a/src/core/actor/DipolarDirectSum.cpp +++ b/src/core/actor/DipolarDirectSum.cpp @@ -31,13 +31,14 @@ if (dipolarDirectSum) { forceActors.remove(dipolarDirectSum); energyActors.remove(dipolarDirectSum); - free(dipolarDirectSum); + delete(dipolarDirectSum); + dipolarDirectSum=NULL; } } -DipolarDirectSum *dipolarDirectSum; +DipolarDirectSum *dipolarDirectSum=0; #endif diff --git a/src/core/actor/DipolarDirectSum.hpp b/src/core/actor/DipolarDirectSum.hpp index e80c2957f6f..eb67448cd08 100644 --- a/src/core/actor/DipolarDirectSum.hpp +++ b/src/core/actor/DipolarDirectSum.hpp @@ -8,6 +8,7 @@ #include "DipolarDirectSum_cuda.hpp" #include "grid.hpp" #include "cuda_interface.hpp" +#include "interaction_data.hpp" #ifndef ACTOR_DIPOLARDIRECTSUM_HPP #define ACTOR_DIPOLARDIRECTSUM_HPP @@ -24,7 +25,10 @@ void DipolarDirectSum_kernel_wrapper_force(dds_float k, int n, float *pos, float class DipolarDirectSum : public Actor { public: DipolarDirectSum(SystemInterface &s) { - if(!s.requestFGpu()) + + k = coulomb.Dbjerrum; + + if(!s.requestFGpu()) std::cerr << "DipolarDirectSum needs access to forces on GPU!" << std::endl; if(!s.requestRGpu()) diff --git a/src/core/actor/DipolarDirectSum_cuda.cu b/src/core/actor/DipolarDirectSum_cuda.cu index 97181192373..3e6d2a515a4 100644 --- a/src/core/actor/DipolarDirectSum_cuda.cu +++ b/src/core/actor/DipolarDirectSum_cuda.cu @@ -136,7 +136,7 @@ __device__ dds_float dipole_ia_energy(int id,dds_float pf, float* r1, float *r2, dds_float dr[3]; dds_float _r1[3],_r2[3],_dip1[3],_dip2[3]; -for (int i=i=0;i<3;i++) +for (int i=0;i<3;i++) { _r1[i]=r1[i]; _r2[i]=r2[i]; @@ -175,7 +175,6 @@ __global__ void DipolarDirectSum_kernel_force(dds_float pf, int i = blockIdx.x * blockDim.x + threadIdx.x; - if(i >= n) return; @@ -215,7 +214,7 @@ __global__ void DipolarDirectSum_kernel_force(dds_float pf, for (int j=i+1;j= n_bonded_ia)) return 3; + // If the bond type to bind particle centers is not a pair bond... + // Check that the bonds have the right number of partners if ((mode & COLLISION_MODE_BOND) && (bonded_ia_params[bond_centers].num != 1)) return 4; + // The bond between the virtual sites can be pair or triple if ((mode & COLLISION_MODE_VS) && !(bonded_ia_params[bond_vs].num == 1 || bonded_ia_params[bond_vs].num == 2)) return 5; + + if (mode & COLLISION_MODE_BIND_THREE_PARTICLES) { + if (bond_three_particles + angle_resolution >= n_bonded_ia) + return 6; + + for (int i = bond_three_particles; i <= bond_three_particles + angle_resolution; i++) { + if (bonded_ia_params[i].num != 2) + return 7; + } + } // Set params collision_params.mode=mode; @@ -89,31 +111,111 @@ int collision_detection_set_params(int mode, double d, int bond_centers, collision_params.bond_vs=bond_vs; collision_params.distance=d; collision_params.vs_particle_type=t; + collision_params.dist_glued_part_to_vs =d2; + collision_params.part_type_to_be_glued =tg; + collision_params.part_type_to_attach_vs_to =tv; + collision_params.part_type_after_glueing =ta; + collision_params.bond_three_particles=bond_three_particles; + collision_params.three_particle_angle_resolution=angle_resolution; + + if (mode & COLLISION_MODE_VS) + make_particle_type_exist(t); + + + if (mode & COLLISION_MODE_GLUE_TO_SURF) + { + make_particle_type_exist(t); + make_particle_type_exist(tg); + make_particle_type_exist(tv); + make_particle_type_exist(ta); + } - make_particle_type_exist(t); mpi_bcast_collision_params(); + + + recalc_forces = 1; + return 0; } +//* Allocate memory for the collision queue / +void prepare_collision_queue() +{ + TRACE(printf("%d: Prepare_collision_queue()\n",this_node)); + number_of_collisions=0; +} + + +bool bond_exists(Particle* p, Particle* partner, int bond_type) +{ + // First check the bonds of p1 + if (p->bl.e) { + int i = 0; + while(i < p->bl.n) { + int size = bonded_ia_params[p->bl.e[i]].num; + + if (p->bl.e[i] == bond_type && + p->bl.e[i + 1] == partner->p.identity) { + // There's a bond, already. Nothing to do for these particles + return true; + } + i += size + 1; + } + } + return false; +} + + +void queue_collision(int part1,int part2, double* point_of_collision) { + + //Get memory for the new entry in the collision queue + number_of_collisions++; + if (number_of_collisions==1) + collision_queue = (collision_struct *) malloc(number_of_collisions*sizeof(collision_struct)); + else + collision_queue = (collision_struct *) realloc (collision_queue,number_of_collisions*sizeof(collision_struct)); + // Save the collision + collision_queue[number_of_collisions-1].pp1 = part1; + collision_queue[number_of_collisions-1].pp2 = part2; + memcpy(collision_queue[number_of_collisions-1].point_of_collision,point_of_collision, 3*sizeof(double)); + + TRACE(printf("%d: Added to queue: Particles %d and %d at %lf %lf %lf\n",this_node,part1,part2,point_of_collision[0],point_of_collision[1],point_of_collision[2])); +} + + // Detect a collision between the given particles. // Add it to the queue in case virtual sites should be added at the point of collision void detect_collision(Particle* p1, Particle* p2) { // The check, whether collision detection is actually turned on is performed in forces.hpp - double dist_betw_part, vec21[3]; int part1, part2, size; + int counts[n_nodes]; + //TRACE(printf("%d: consider particles %d and %d\n", this_node, p1->p.identity, p2->p.identity)); - TRACE(printf("%d: consider particles %d and %d\n", this_node, p1->p.identity, p2->p.identity)); - + double vec21[3]; // Obtain distance between particles - dist_betw_part = sqrt(distance2vec(p1->r.p, p2->r.p, vec21)); + double dist_betw_part = sqrt(distance2vec(p1->r.p, p2->r.p, vec21)); + //TRACE(printf("%d: Distance between particles %lf %lf %lf, Scalar: %f\n",this_node,vec21[0],vec21[1],vec21[2], dist_betw_part)); if (dist_betw_part > collision_params.distance) return; - TRACE(printf("%d: particles %d and %d on bonding distance %lf\n", this_node, p1->p.identity, p2->p.identity, dist_betw_part)); + //TRACE(printf("%d: particles %d and %d within bonding distance %lf\n", this_node, p1->p.identity, p2->p.identity, dist_betw_part)); + // If we are in the glue to surface mode, check that the particles + // are of the right type + if (collision_params.mode & COLLISION_MODE_GLUE_TO_SURF) { + if (! ( + ((p1->p.type==collision_params.part_type_to_be_glued) + && (p2->p.type ==collision_params.part_type_to_attach_vs_to)) + || + ((p2->p.type==collision_params.part_type_to_be_glued) + && (p1->p.type ==collision_params.part_type_to_attach_vs_to))) + ) { + return; + } + } part1 = p1->p.identity; part2 = p2->p.identity; @@ -129,98 +231,180 @@ void detect_collision(Particle* p1, Particle* p2) return; #endif - // Check, if there's already a bond between the particles - // First check the bonds of p1 - if (p1->bl.e) { - int i = 0; - while(i < p1->bl.n) { - size = bonded_ia_params[p1->bl.e[i]].num; - - if (p1->bl.e[i] == collision_params.bond_centers && - p1->bl.e[i + 1] == part2) { - // There's a bond, already. Nothing to do for these particles - return; - } - i += size + 1; - } - } - if (p2->bl.e) { - // Check, if a bond is already stored in p2 - int i = 0; - while(i < p2->bl.n) { - size = bonded_ia_params[p2->bl.e[i]].num; + if (p1==p2) + return; - /* COMPARE P2 WITH P1'S BONDED PARTICLES*/ + // Check, if there's already a bond between the particles + if (bond_exists(p1,p2, collision_params.bond_centers)) + return; + + if (bond_exists(p2,p1, collision_params.bond_centers)) + return; - if (p2->bl.e[i] == collision_params.bond_centers && - p2->bl.e[i + 1] == part1) { - return; - } - i += size + 1; - } - } TRACE(printf("%d: no previous bond, binding\n", this_node)); /* If we're still here, there is no previous bond between the particles, we have a new collision */ - /* create marking bond between the colliding particles immediately */ if (collision_params.mode & COLLISION_MODE_BOND) { - int bondG[2]; - int primary = part1, secondary = part2; - // put the bond to the physical particle; at least one partner always is - if (p1->l.ghost) { - primary = part2; - secondary = part1; + + // do not create bond between ghost particles + if (p1->l.ghost && p2->l.ghost) { + TRACE(printf("Both particles %d and %d are ghost particles", p1->p.identity, p2->p.identity)); + return; } - bondG[0]=collision_params.bond_centers; - bondG[1]=secondary; - local_change_bond(primary, bondG, 0); - } - if (collision_params.mode & (COLLISION_MODE_VS | COLLISION_MODE_EXCEPTION)) { - /* If we also create virtual sites or throw an exception, we add the collision + + double new_position[3]; + /* If we also create virtual sites or bind three particles, or throw an exception, we add the collision to the queue to process later */ - // Point of collision - double new_position[3]; - for (int i=0;i<3;i++) { - new_position[i] = p1->r.p[i] - vec21[i] * 0.50; + double c; + + // If not in the glue_to_surface-mode, the point of collision + // is in the middle of the vecotr connecting the particle + // centers + if (! (collision_params.mode & COLLISION_MODE_GLUE_TO_SURF)) + c=0.5; + else + { + // Find out, which is the particle to be glued. + // Swap particle, if need be + if ((p1->p.type==collision_params.part_type_to_be_glued) + && (p2->p.type ==collision_params.part_type_to_attach_vs_to)) + { + c = collision_params.dist_glued_part_to_vs/dist_betw_part; + } + else if ((p2->p.type==collision_params.part_type_to_be_glued) + && (p1->p.type ==collision_params.part_type_to_attach_vs_to)) + { + // we swap the particle ids so that the virtual site is always attached to p2 + int tmp=part1; + part1=part2; + part2=tmp; + // We need the negative sign for the prefactor, as we di + // we did not flip the vec21 when swapping particles + c = -collision_params.dist_glued_part_to_vs/dist_betw_part; + } + else + { + printf("Something is wrong %s %d\n",__FILE__,":"+__LINE__); + } + } + for (int i=0;i<3;i++) { + new_position[i] = p1->r.p[i] - vec21[i] * c; } - - number_of_collisions++; - // Allocate mem for the new collision info - collision_queue = (collision_struct *) Utils::realloc (collision_queue, (number_of_collisions) * sizeof(collision_struct)); - - // Save the collision - collision_queue[number_of_collisions-1].pp1 = part1; - collision_queue[number_of_collisions-1].pp2 = part2; - for (int i=0;i<3;i++) { - collision_queue[number_of_collisions-1].point_of_collision[i] = new_position[i]; - } + + + queue_collision(part1,part2,new_position); } } -void prepare_collision_queue() + + +// Considers three particles for three_particle_binding and performs +// the binding if the criteria are met // +void coldet_do_three_particle_bond(Particle* p, Particle* p1, Particle* p2) { - - number_of_collisions=0; + double vec21[3]; + // If p1 and p2 are not closer or equal to the cutoff distance, skip + // p1: + get_mi_vector(vec21,p->r.p,p1->r.p); + if (sqrt(sqrlen(vec21)) > collision_params.distance) + return; + // p2: + get_mi_vector(vec21,p->r.p,p2->r.p); + if (sqrt(sqrlen(vec21)) > collision_params.distance) + return; + + //TRACE(printf("%d: checking three particle bond %d %d %d\n", this_node, p1->p.identity, p->p.identity, p2->p.identity)); + + // Check, if there already is a three-particle bond centered on p + // with p1 and p2 as partners. If so, skip this triplet. + // Note that the bond partners can appear in any order. + + // Iterate over existing bonds of p + + if (p->bl.e) { + int b = 0; + while (b < p->bl.n) { + int size = bonded_ia_params[p->bl.e[b]].num; + + //TRACE(printf("%d:--1-- checking bond of type %d and length %d of particle %d\n", this_node, p->bl.e[b], bonded_ia_params[p->bl.e[b]].num, p->p.identity)); + + if (size==2) { + // Check if the bond type is within the range used by the collision detection, + if ((p->bl.e[b] >= collision_params.bond_three_particles) & (p->bl.e[b] <=collision_params.bond_three_particles + collision_params.three_particle_angle_resolution)) { + // check, if p1 and p2 are the bond partners, (in any order) + // if yes, skip triplet + if ( + ((p->bl.e[b+1]==p1->p.identity) && (p->bl.e[b+2] ==p2->p.identity)) + || + ((p->bl.e[b+1]==p2->p.identity) && (p->bl.e[b+2] ==p1->p.identity)) + ) + return; + } // if bond type + } // if size==2 + + // Go to next bond + b += size + 1; + } // bond loop + } // if bond list defined - collision_queue = (collision_struct *) Utils::malloc (sizeof(collision_struct)); + //TRACE(printf("%d: proceeding to install three particle bond %d %d %d\n", this_node, p1->p.identity, p->p.identity, p2->p.identity)); + // If we are still here, we need to create angular bond + // First, find the angle between the particle p, p1 and p2 + double cosine=0.0; + + double vec1[3],vec2[3]; + /* vector from p to p1 */ + get_mi_vector(vec1, p->r.p, p1->r.p); + // Normalize + double dist2 = sqrlen(vec1); + double d1i = 1.0 / sqrt(dist2); + for(int j=0;j<3;j++) vec1[j] *= d1i; + + /* vector from p to p2 */ + get_mi_vector(vec2, p->r.p, p2->r.p); + // normalize + dist2 = sqrlen(vec2); + double d2i = 1.0 / sqrt(dist2); + for(int j=0;j<3;j++) vec2[j] *= d2i; + + /* scalar produvt of vec1 and vec2 */ + cosine = scalar(vec1, vec2); + + // Handle case where cosine is nearly 1 or nearly -1 + if ( cosine > TINY_COS_VALUE) + cosine = TINY_COS_VALUE; + if ( cosine < -TINY_COS_VALUE) + cosine = -TINY_COS_VALUE; + + // Bond angle + double phi = acos(cosine); + + // We find the bond id by dividing the range from 0 to pi in + // three_particle_angle_resolution steps and by adding the id + // of the bond for zero degrees. + int bond_id =floor(phi/M_PI * (collision_params.three_particle_angle_resolution-1) +0.5) +collision_params.bond_three_particles; + + // Create the bond + + // First, fill bond data structure + int bondT[3]; + bondT[0] = bond_id; + bondT[1] = p1->p.identity; + bondT[2] = p2->p.identity; + local_change_bond(p->p.identity, bondT, 0); } -// Handle the collisions stored in the queue -void handle_collisions () +// If activated, throws an exception for each collision which can be +// parsed by the script interface +void handle_exception_throwing_for_single_collision(int i) { - // printf("number of collisions in handle collision are %d\n",number_of_collisions); - - for (int i = 0; i < number_of_collisions; i++) { - // printf("Handling collision of particles %d %d\n", collision_queue[i].pp1, collision_queue[i].pp2); - // fflush(stdout); - if (collision_params.mode & (COLLISION_MODE_EXCEPTION)) { int id1, id2; @@ -236,57 +420,342 @@ void handle_collisions () msg << "collision between particles " << id1 << " and " <p.isVirtual=1; + #ifdef ROTATION_PER_PARTICLE + (local_particles[relate_to])->p.rotation=14; + #endif + (local_particles[max_seen_particle])->p.type=collision_params.vs_particle_type; +} + + +void bind_at_poc_create_bond_between_vs(int i) +{ + int bondG[3]; + + switch (bonded_ia_params[collision_params.bond_vs].num) { + case 1: { + // Create bond between the virtual particles + bondG[0] = collision_params.bond_vs; + bondG[1] = max_seen_particle-1; + local_change_bond(max_seen_particle, bondG, 0); + break; + } + case 2: { + // Create 1st bond between the virtual particles + bondG[0] = collision_params.bond_vs; + bondG[1] = collision_queue[i].pp1; + bondG[2] = collision_queue[i].pp2; + local_change_bond(max_seen_particle, bondG, 0); + local_change_bond(max_seen_particle-1, bondG, 0); + break; + } + } +} + +void glue_to_surface_bind_vs_to_pp1(int i) +{ + int bondG[3]; + // Create bond between the virtual particles + bondG[0] = collision_params.bond_vs; + bondG[1] = max_seen_particle; + local_change_bond(collision_queue[i].pp1, bondG, 0); + local_particles[collision_queue[i].pp1]->p.type=collision_params.part_type_after_glueing; +} + +#endif + +void gather_collision_queue(int* counts) +{ + int displacements[n_nodes]; // offsets into collisions + + // Initialize number of collisions gathered from all processors + for (int a=0;an; a++) { + p=&cell->part[a]; + // for all p: + for (int ij=0; ijp.isVirtual=1; - (local_particles[max_seen_particle])->p.type=collision_params.vs_particle_type; + // Check, whether p is equal to one of the particles in the + // collision. If so, skip + if ((p->p.identity ==p1->p.identity) || ( p->p.identity == p2->p.identity)) { + continue; + } - // Virtual particle related to 2nd particle of the collision - place_particle(max_seen_particle+1,collision_queue[i].point_of_collision); - vs_relate_to(max_seen_particle,collision_queue[i].pp2); - (local_particles[max_seen_particle])->p.isVirtual=1; - (local_particles[max_seen_particle])->p.type=collision_params.vs_particle_type; + // The following checks, + // if the particle p is closer that the cutoff from p1 and/or p2. + // If yes, three particle bonds are created on all particles + // which have two other particles within the cutoff distance, + // unless such a bond already exists + + // We need all cyclical permutations, here + // (bond is placed on 1st particle, order of bond partners + // does not matter, so we don't neet non-cyclic permutations): + coldet_do_three_particle_bond(p,p1,p2); + coldet_do_three_particle_bond(p1,p,p2); + coldet_do_three_particle_bond(p2,p,p1); - switch (bonded_ia_params[collision_params.bond_vs].num) { - case 1: { - // Create bond between the virtual particles - bondG[0] = collision_params.bond_vs; - bondG[1] = max_seen_particle-1; - local_change_bond(max_seen_particle, bondG, 0); - break; - } - case 2: { - // Create 1st bond between the virtual particles - bondG[0] = collision_params.bond_vs; - bondG[1] = collision_queue[i].pp1; - bondG[2] = collision_queue[i].pp2; - local_change_bond(max_seen_particle, bondG, 0); - local_change_bond(max_seen_particle-1, bondG, 0); - break; - } + } + } + } +} + + +// Goes through the collision queue and for each pair in it +// looks for a third particle by using the domain decomposition +// cell system. If found, it performs three particle binding +void three_particle_binding_domain_decomposition() +{ + // We have domain decomposition + + // Indices of the cells in which the colliding particles reside + int cellIdx[2][3]; + + // Iterate over collision queue + + for (int id=0;idr.p,cellIdx[0]); + dd_position_to_cell_indices(p2->r.p,cellIdx[1]); + + // Iterate over the cells + their neighbors + // if p1 and p2 are in the same cell, we don't need to consider it 2x + int lim=1; + + if ((cellIdx[0][0]==cellIdx[1][0]) && (cellIdx[0][1]==cellIdx[1][1]) && (cellIdx[0][2]==cellIdx[1][2])) + lim=0; // Only consider the 1st cell + + for (int j=0;j<=lim;j++) { + + // Iterate the cell with indices cellIdx[j][] and all its neighbors. + // code taken from dd_init_cell_interactions() + for(int p=cellIdx[j][0]-1; p<=cellIdx[j][0]+1; p++) + for(int q=cellIdx[j][1]-1; q<=cellIdx[j][1]+1; q++) + for(int r=cellIdx[j][2]-1; r<=cellIdx[j][2]+1; r++) { + int ind2 = get_linear_index(p,q,r,dd.ghost_cell_grid); + Cell* cell=cells+ind2; + + // Iterate over particles in this cell + for(int a=0; an; a++) { + Particle* P=&cell->part[a]; + // for all p: + // Check, whether p is equal to one of the particles in the + // collision. If so, skip + if ((P->p.identity ==p1->p.identity) || (P->p.identity == p2->p.identity)) { + //TRACE(printf("same particle\n")); + continue; + } + + // The following checks, + // if the particle p is closer that the cutoff from p1 and/or p2. + // If yes, three particle bonds are created on all particles + // which have two other particles within the cutoff distance, + // unless such a bond already exists + + // We need all cyclical permutations, here + // (bond is placed on 1st particle, order of bond partners + // does not matter, so we don't need non-cyclic permutations): + + if (P->l.ghost) { + //TRACE(printf("%d: center particle is ghost: %d\n", this_node, P->p.identity)); + continue; + } + //TRACE(printf("%d: LOOP: %d Handling collision of particles FIRST CONFIGURATION %d %d %d\n", this_node, id, p1->p.identity, P->p.identity, p2->p.identity)); + coldet_do_three_particle_bond(P,p1,p2); + + if (p1->l.ghost) { + //TRACE(printf("%d: center particle is ghost: %d\n", this_node, p1->p.identity)); + continue; + } + + coldet_do_three_particle_bond(p1,P,p2); + + if (p2->l.ghost) { + //TRACE(printf("%d: center particle is ghost: %d\n", this_node, p2->p.identity)); + continue; + } + + coldet_do_three_particle_bond(p2,P,p1); + + } // loop over particles in this cell + + } // Loop over cell + + } // Loop over particles if they are in different cells + + } // If local particles exist + + } // Loop over total collisions +} + + +// Handle the collisions stored in the queue +void handle_collisions () +{ + + TRACE(printf("%d: handle_collisions: number of collisions in queue %d\n",this_node,number_of_collisions)); + + if (collision_params.mode & COLLISION_MODE_EXCEPTION) + for (int i=0;il.ghost) { + primary = collision_queue[i].pp2; + secondary = collision_queue[i].pp1; + TRACE(printf("%d: particle-%d is ghost", this_node, collision_queue[i].pp1)); } + int bondG[2]; + bondG[0]=collision_params.bond_centers; + bondG[1]=secondary; + local_change_bond(primary, bondG, 0); + TRACE(printf("%d: Adding bond %d->%d\n",this_node, primary,secondary)); } -#endif } +#ifdef VIRTUAL_SITES_RELATIVE + // If one of the collision modes is active which places virtual sites, we go over the queue to handle them + if ((collision_params.mode & COLLISION_MODE_VS) || (collision_params.mode & COLLISION_MODE_GLUE_TO_SURF)) { + for (int i=0;i0) { + + // If we don't have domain decomposition, we need to do a full sweep over all + // particles in the system. (slow) + if (cell_structure.type!=CELL_STRUCTURE_DOMDEC) { + three_particle_binding_full_search(); + } // if cell structure != domain decomposition + else + { + three_particle_binding_domain_decomposition(); + } // If we have doamin decomposition + + } // if number of collisions of this node > 0 + + if (total_collisions>0) + free(gathered_queue); + total_collisions = 0; + } // if TPB + + // If a collision method is active which places particles, resorting might be needed + TRACE(printf("%d: Resort particles is %d\n",this_node,resort_particles)); + if (collision_params.mode & (COLLISION_MODE_VS | COLLISION_MODE_GLUE_TO_SURF)) + { + // NOTE!! this has to be changed to total_collisions, once parallelization + // is implemented + + if (number_of_collisions >0) + { + announce_resort_particles(); + } + } + // Reset the collision queue - if(collision_queue && (number_of_collisions > 0)) { + if (number_of_collisions>0) free(collision_queue); - collision_queue = 0; - number_of_collisions = 0; - } + + number_of_collisions = 0; + - announce_resort_particles(); } #endif diff --git a/src/core/collision.hpp b/src/core/collision.hpp index d55d49422ac..07c6886c91c 100755 --- a/src/core/collision.hpp +++ b/src/core/collision.hpp @@ -44,6 +44,10 @@ together. This prevents the particles from sliding against each other. Requires VIRTUAL_SITES_RELATIVE and COLLISION_MODE_BOND*/ #define COLLISION_MODE_VS 4 +/** Glue a particle to a speciffic spot on the surface of an other */ +#define COLLISION_MODE_GLUE_TO_SURF 8 +/// Three particle binding mode +#define COLLISION_MODE_BIND_THREE_PARTICLES 16 /*@}*/ typedef struct { @@ -57,6 +61,20 @@ typedef struct { int mode; /// distance at which particles are bound double distance; + /// For mode "glue to surface": The distance from the particle which is to be glued to the new virtual site + double dist_glued_part_to_vs; + /// For mode "glue to surface": The particle type being glued + int part_type_to_be_glued; + /// For mode "glue to surface": The particle type to which the virtual site is attached + int part_type_to_attach_vs_to; + /// Particle type to which the newly glued particle is converd + int part_type_after_glueing; + /// first bond type (for zero degrees)) used for the three-particle bond (angle potential) + int bond_three_particles; + /// Number of angle bonds to use (angular resolution) + /// different angle bonds with different equilibrium angles + /// Are expected to have ids immediately following to bond_three_particles + int three_particle_angle_resolution; } Collision_parameters; /// Parameters for collision detection @@ -80,8 +98,13 @@ void handle_collisions(); @param bond_vs is the type of the bond between the virtual particles, if using noslip bonds @param t is the type of the virtual sites, if using noslip bonds + @param d2 for the "glue to surface" mode is the distance between the particle to be glued and the new virtual site + @param tg for the "glue to surface" is the type of the particle being glued + @param tv for the "glue to surface" is the type of the particle to which the virtual site is attached + @param bond_three_particles is the three-particle-bond parameter + @param angle_resolution is the three_particle_angle_resolution parameter in order to define different angle bonds */ -int collision_detection_set_params(int mode, double d, int bond_centers, int bond_vs,int t); +int collision_detection_set_params(int mode, double d, int bond_centers, int bond_vs,int t,int d2, int tg, int tv, int ta, int bond_three_particles, int angle_resolution); #endif diff --git a/src/core/domain_decomposition.cpp b/src/core/domain_decomposition.cpp index 97e1df6ad44..87da6417cb0 100755 --- a/src/core/domain_decomposition.cpp +++ b/src/core/domain_decomposition.cpp @@ -646,6 +646,25 @@ Cell *dd_position_to_cell(double pos[3]) return &cells[i]; } +void dd_position_to_cell_indices(double pos[3],int* idx) +{ + int i; + double lpos; + + for(i=0;i<3;i++) { + lpos = pos[i] - my_left[i]; + + idx[i] = (int)(lpos*dd.inv_cell_size[i])+1; + + if (idx[i] < 1) { + idx[i] = 1; + } + else if (idx[i] > dd.cell_grid[i]) { + idx[i] = dd.cell_grid[i]; + } + } +} + /*************************************************/ /** Append the particles in pl to \ref local_cells and update \ref local_particles. diff --git a/src/core/domain_decomposition.hpp b/src/core/domain_decomposition.hpp index b9ca24ee367..23e32a17f5c 100755 --- a/src/core/domain_decomposition.hpp +++ b/src/core/domain_decomposition.hpp @@ -198,6 +198,9 @@ void dd_exchange_and_sort_particles(int global_flag); /** implements \ref CellStructure::position_to_cell. */ Cell *dd_position_to_cell(double pos[3]); +/** Get three cell indices (coordinates in cell gird) from particle position */ +void dd_position_to_cell_indices(double pos[3],int* idx); + /** calculate physical (processor) minimal number of cells */ int calc_processor_min_num_cells(); diff --git a/src/core/forces_inline.hpp b/src/core/forces_inline.hpp index d262e2b73bb..2d66425afae 100644 --- a/src/core/forces_inline.hpp +++ b/src/core/forces_inline.hpp @@ -320,9 +320,6 @@ espressoSystemInterface.update(); // mark that forces are now up-to-date recalc_forces = 0; -#ifdef COLLISION_DETECTION - handle_collisions(); -#endif } diff --git a/src/core/ghosts.cpp b/src/core/ghosts.cpp index 12017d8447f..f7bb191ca3c 100755 --- a/src/core/ghosts.cpp +++ b/src/core/ghosts.cpp @@ -448,10 +448,10 @@ void cell_cell_transfer(GhostCommunication *gc, int data_parts) memmove(&pt2->p, &pt1->p, sizeof(ParticleProperties)); #ifdef GHOSTS_HAVE_BONDS realloc_intlist(&(pt2->bl), pt2->bl.n = pt1->bl.n); - memmove(&pt2->bl.e, &pt1->bl.e, pt1->bl.n*sizeof(int)); + memmove(pt2->bl.e, pt1->bl.e, pt1->bl.n*sizeof(int)); #ifdef EXCLUSIONS realloc_intlist(&(pt2->el), pt2->el.n = pt1->el.n); - memmove(&pt2->el.e, &pt1->el.e, pt1->el.n*sizeof(int)); + memmove(pt2->el.e, pt1->el.e, pt1->el.n*sizeof(int)); #endif #endif } @@ -718,6 +718,7 @@ void invalidate_ghosts() particle array. */ if( &(part[p]) == local_particles[part[p].p.identity] ) local_particles[part[p].p.identity] = NULL; + free_particle(part+p); } ghost_cells.cell[c]->n = 0; } diff --git a/src/core/integrate.cpp b/src/core/integrate.cpp index 5e4245b784c..dddfcb37c1b 100755 --- a/src/core/integrate.cpp +++ b/src/core/integrate.cpp @@ -553,6 +553,9 @@ void integrate_vv(int n_steps, int reuse_forces) /* Propagate time: t = t+dt */ sim_time += time_step; } + #ifdef COLLISION_DETECTION + handle_collisions(); + #endif } /* verlet list statistics */ diff --git a/src/core/mpifake/mpi.h b/src/core/mpifake/mpi.h index 545803f8596..5d0e4309075 100755 --- a/src/core/mpifake/mpi.h +++ b/src/core/mpifake/mpi.h @@ -195,6 +195,13 @@ inline int MPI_Allgather(void *sbuf, int scount, MPI_Datatype sdtype, void *rbuf, int rcount, MPI_Datatype rdtype, MPI_Comm comm) { return mpifake_sendrecv(sbuf, scount, sdtype, rbuf, rcount, rdtype); } +inline int MPI_Allgatherv(void *sbuf, int scount, MPI_Datatype sdtype, + void *rbuf, int* rcount, int* displ, MPI_Datatype rdtype, + MPI_Comm comm) +{ return mpifake_sendrecv(sbuf, scount, sdtype, rbuf, rcount[0], rdtype); } + + + inline int MPI_Scatter(void *sbuf, int scount, MPI_Datatype sdtype, void *rbuf, int rcount, MPI_Datatype rdtype, int root, MPI_Comm comm) diff --git a/src/core/virtual_sites_relative.cpp b/src/core/virtual_sites_relative.cpp index 1518c298249..515f4fb3430 100755 --- a/src/core/virtual_sites_relative.cpp +++ b/src/core/virtual_sites_relative.cpp @@ -284,6 +284,8 @@ int vs_relate_to(int part_num, int relate_to) quat[0]=1; quat[1]=quat[2]=quat[3]=0; } + free_particle(&p_relate_to); + free_particle(&p_current); // Set the particle id of the particle we want to relate to, the distnace // and the relative orientation diff --git a/src/tcl/collision_tcl.cpp b/src/tcl/collision_tcl.cpp index 1fa1ac78d59..a5c7178309e 100755 --- a/src/tcl/collision_tcl.cpp +++ b/src/tcl/collision_tcl.cpp @@ -43,19 +43,39 @@ int tclcommand_on_collision(ClientData data, Tcl_Interp *interp, int argc, char /* this one can be combined with the rest */ if (collision_params.mode & COLLISION_MODE_EXCEPTION) { sprintf(s, " exception"); + Tcl_AppendResult(interp, s + 1, (char*) NULL); } if (collision_params.mode & COLLISION_MODE_VS) { sprintf(s, " bind_at_point_of_collision %f %d %d %d", collision_params.distance, collision_params.bond_centers, collision_params.bond_vs, collision_params.vs_particle_type); + Tcl_AppendResult(interp, s + 1, (char*) NULL); + return TCL_OK; + } + + if (collision_params.mode & COLLISION_MODE_BIND_THREE_PARTICLES) { + sprintf(s, " bind_three_particles %f %d %d %d", + collision_params.distance, collision_params.bond_centers, + collision_params.bond_three_particles, collision_params.three_particle_angle_resolution); + Tcl_AppendResult(interp, s + 1, (char*) NULL); + } + else if (collision_params.mode & COLLISION_MODE_GLUE_TO_SURF) { + sprintf(s, " glue_to_surface %f %d %d %d %d %d %d %f", + collision_params.distance, collision_params.bond_centers, + collision_params.bond_vs, collision_params.vs_particle_type, + collision_params.part_type_to_be_glued, + collision_params.part_type_to_attach_vs_to, + collision_params.part_type_after_glueing, + collision_params.dist_glued_part_to_vs); + Tcl_AppendResult(interp, s + 1, (char*) NULL); } else if (collision_params.mode & COLLISION_MODE_BOND) { sprintf(s, " bind_centers %f %d", collision_params.distance, collision_params.bond_centers); + Tcl_AppendResult(interp, s + 1, (char*) NULL); } // first character is always the separating space - Tcl_AppendResult(interp, s + 1, (char*) NULL); return TCL_OK; } @@ -63,16 +83,26 @@ int tclcommand_on_collision(ClientData data, Tcl_Interp *interp, int argc, char // Otherwise, we set parameters if (ARG0_IS_S("off")) { - collision_detection_set_params(0,0,0,0,0); + collision_detection_set_params(0,0,0,0,0,0,0,0,0,0,0); return TCL_OK; } else { /* parameters of collision_detection_set_params */ int mode = 0; - double d = 0; + + // Distances + double d,d2 = 0; + + // Bond types int bond_centers = 0; int bond_vs = 0; - int t = 0; + + // Particle types for virtual sites based based methods + int t,tg,tv,ta = 0; + + // /bond types for three particle binding + int bond_three_particles=0; + int angle_resolution=0; if (ARG0_IS_S("exception")) { mode = COLLISION_MODE_EXCEPTION; @@ -122,12 +152,76 @@ int tclcommand_on_collision(ClientData data, Tcl_Interp *interp, int argc, char } argc -= 5; argv += 5; } + else if (ARG0_IS_S("glue_to_surface")) { + mode |= COLLISION_MODE_BOND | COLLISION_MODE_GLUE_TO_SURF; + if (argc != 9) { + Tcl_AppendResult(interp, "Not enough parameters, need a distance, two bond types, four particle types and another distance as args.", (char*) NULL); + return TCL_ERROR; + } + if (!ARG_IS_D(1,d)) { + Tcl_AppendResult(interp, "Need a distance as 1st arg.", (char*) NULL); + return TCL_ERROR; + } + if (!ARG_IS_I(2,bond_centers)) { + Tcl_AppendResult(interp, "Need a bond type as 2nd arg.", (char*) NULL); + return TCL_ERROR; + } + if (!ARG_IS_I(3,bond_vs)) { + Tcl_AppendResult(interp, "Need a bond type as 3rd arg.", (char*) NULL); + return TCL_ERROR; + } + if (!ARG_IS_I(4,t)) { + Tcl_AppendResult(interp, "Need a particle type as 4th arg.", (char*) NULL); + return TCL_ERROR; + } + if (!ARG_IS_I(5,tg)) { + Tcl_AppendResult(interp, "Need a particle type as 5th arg.", (char*) NULL); + return TCL_ERROR; + } + if (!ARG_IS_I(6,tv)) { + Tcl_AppendResult(interp, "Need a particle type as 6th arg.", (char*) NULL); + return TCL_ERROR; + } + if (!ARG_IS_I(7,ta)) { + Tcl_AppendResult(interp, "Need a particle type as 7th arg.", (char*) NULL); + return TCL_ERROR; + } + if (!ARG_IS_D(8,d2)) { + Tcl_AppendResult(interp, "Need a distance as 8th arg.", (char*) NULL); + return TCL_ERROR; + } + argc -= 9; argv += 8; + } + else if (ARG0_IS_S("bind_three_particles")) { + mode |= COLLISION_MODE_BIND_THREE_PARTICLES | COLLISION_MODE_BOND; + if (argc != 5) { + Tcl_AppendResult(interp, "Not enough parameters, need a distance and two bond types.", (char*) NULL); + return TCL_ERROR; + } + if (!ARG_IS_D(1,d)) { + Tcl_AppendResult(interp, "Need a distance as 1st arg.", (char*) NULL); + return TCL_ERROR; + } + if (!ARG_IS_I(2,bond_centers)) { + Tcl_AppendResult(interp, "Need a bond type as 2nd arg.", (char*) NULL); + return TCL_ERROR; + } + if (!ARG_IS_I(3,bond_three_particles)) { + Tcl_AppendResult(interp, "Need a bond type as 3rd arg.", (char*) NULL); + return TCL_ERROR; + } + if (!ARG_IS_I(4,angle_resolution)) { + Tcl_AppendResult(interp, "Need an angle resolution as 4th arg.", (char*) NULL); + return TCL_ERROR; + } + argc -= 5; argv += 5; + } else { Tcl_AppendResult(interp, "\"", argv[0], "\" is not a valid collision detection mode.", (char*) NULL); return TCL_ERROR; } - int res = collision_detection_set_params(mode,d,bond_centers,bond_vs,t); + int res = collision_detection_set_params(mode,d,bond_centers,bond_vs,t,d2,tg,tv,ta,bond_three_particles,angle_resolution); switch (res) { case 1: @@ -145,6 +239,12 @@ int tclcommand_on_collision(ClientData data, Tcl_Interp *interp, int argc, char case 5: Tcl_AppendResult(interp, "Virtual particles need a pair bond or triple bond.", (char*) NULL); return TCL_ERROR; + case 6: + Tcl_AppendResult(interp, "Not enough angular bonds.", (char*) NULL); + return TCL_ERROR; + case 7: + Tcl_AppendResult(interp, "bond_three_particles needs triple bonds.", (char*) NULL); + return TCL_ERROR; } return TCL_OK; diff --git a/src/tcl/interaction_data_tcl.cpp b/src/tcl/interaction_data_tcl.cpp index 76149c91ee2..17f0f57244a 100755 --- a/src/tcl/interaction_data_tcl.cpp +++ b/src/tcl/interaction_data_tcl.cpp @@ -1090,7 +1090,6 @@ int tclcommand_inter_parse_rest(Tcl_Interp * interp, int argc, char ** argv) #endif } - Tcl_AppendResult(interp, "unknown interaction type \"", argv[0], "\"", (char *) NULL); diff --git a/src/tcl/statistics_tcl.cpp b/src/tcl/statistics_tcl.cpp index 1de09492d87..e404c7503ae 100755 --- a/src/tcl/statistics_tcl.cpp +++ b/src/tcl/statistics_tcl.cpp @@ -2737,39 +2737,48 @@ int tclcommand_analyze(ClientData data, Tcl_Interp *interp, int argc, char **arg else if (ARG1_IS_S(name)) err = parser(interp, argc - 2, argv + 2) #define REGISTER_ANALYSIS_W_ARG(name, parser, arg) \ else if (ARG1_IS_S(name)) err = parser(interp, arg, argc - 2, argv + 2) - +#define REGISTER_ANALYSIS_WARN(name, parser) \ + else if (ARG1_IS_S(name)) { \ + fprintf(stderr, "WARNING: analyze %s may be deprecated in the future if no documentation or test cases are provided. See open issues at: github.com/espressomd/espresso/issues/ .\n", name); \ + err = parser(interp, argc - 2, argv + 2); \ + } +#define REGISTER_ANALYSIS_W_ARG_WARN(name, parser, arg) \ + else if (ARG1_IS_S(name)) { \ + fprintf(stderr, "WARNING: analyze %s may be deprecated in the future if no documentation or test cases are provided. See open issues at: github.com/espressomd/espresso/issues/ .\n", name); \ + err = parser(interp, arg, argc - 2, argv + 2); \ + } /* for the elses below */ if (0); REGISTER_ANALYZE_OPTION("set", tclcommand_analyze_parse_set); #if defined(LB) || defined(LB_GPU) REGISTER_ANALYZE_OPTION("fluid", tclcommand_analyze_parse_fluid); #endif - REGISTER_ANALYSIS("get_folded_positions", tclcommand_analyze_parse_get_folded_positions); - REGISTER_ANALYSIS("wallstuff", tclcommand_analyze_wallstuff); - REGISTER_ANALYSIS("current", tclcommand_analyze_current); + REGISTER_ANALYSIS_WARN("get_folded_positions", tclcommand_analyze_parse_get_folded_positions) + REGISTER_ANALYSIS_WARN("wallstuff", tclcommand_analyze_wallstuff) + REGISTER_ANALYSIS_WARN("current", tclcommand_analyze_current) #ifdef MODES REGISTER_ANALYZE_OPTION("set_bilayer", tclcommand_analyze_parse_bilayer_set); - REGISTER_ANALYSIS("modes2d", tclcommand_analyze_parse_modes2d); - REGISTER_ANALYSIS("bilayer_density_profile", tclcommand_analyze_parse_bilayer_density_profile); + REGISTER_ANALYSIS_WARN("modes2d", tclcommand_analyze_parse_modes2d) + REGISTER_ANALYSIS_WARN("bilayer_density_profile", tclcommand_analyze_parse_bilayer_density_profile) REGISTER_ANALYSIS("cylindrical_average", tclcommand_analyze_parse_cylindrical_average); - REGISTER_ANALYSIS("radial_density_map", tclcommand_analyze_parse_radial_density_map); - REGISTER_ANALYSIS("get_lipid_orients", tclcommand_analyze_parse_get_lipid_orients); - REGISTER_ANALYSIS("lipid_orient_order", tclcommand_analyze_parse_lipid_orient_order); + REGISTER_ANALYSIS_WARN("radial_density_map", tclcommand_analyze_parse_radial_density_map) + REGISTER_ANALYSIS_WARN("get_lipid_orients", tclcommand_analyze_parse_get_lipid_orients) + REGISTER_ANALYSIS_WARN("lipid_orient_order", tclcommand_analyze_parse_lipid_orient_order) #endif - REGISTER_ANALYSIS("mol", tclcommand_analyze_parse_mol); - REGISTER_ANALYSIS("cluster_size_dist", tclcommand_analyze_parse_cluster_size_dist); + REGISTER_ANALYSIS_WARN("mol", tclcommand_analyze_parse_mol) + REGISTER_ANALYSIS_WARN("cluster_size_dist", tclcommand_analyze_parse_cluster_size_dist) REGISTER_ANALYSIS("mindist", tclcommand_analyze_parse_mindist); - REGISTER_ANALYSIS("aggregation", tclcommand_analyze_parse_aggregation); + REGISTER_ANALYSIS_WARN("aggregation", tclcommand_analyze_parse_aggregation) REGISTER_ANALYSIS("centermass", tclcommand_analyze_parse_centermass); - REGISTER_ANALYSIS("angularmomentum", tclcommand_analyze_parse_angularmomentum); - REGISTER_ANALYSIS("MSD", tclcommand_analyze_parse_MSD); - REGISTER_ANALYSIS("dipmom_normal", tclcommand_analyze_parse_and_print_dipole); - REGISTER_ANALYSIS("momentofinertiamatrix", tclcommand_analyze_parse_momentofinertiamatrix); - REGISTER_ANALYSIS("gyration_tensor", tclcommand_analyze_parse_gyration_tensor); - REGISTER_ANALYSIS("find_principal_axis", tclcommand_analyze_parse_find_principal_axis); + REGISTER_ANALYSIS_WARN("angularmomentum", tclcommand_analyze_parse_angularmomentum) + REGISTER_ANALYSIS_WARN("MSD", tclcommand_analyze_parse_MSD) + REGISTER_ANALYSIS_WARN("dipmom_normal", tclcommand_analyze_parse_and_print_dipole) + REGISTER_ANALYSIS_WARN("momentofinertiamatrix", tclcommand_analyze_parse_momentofinertiamatrix) + REGISTER_ANALYSIS_WARN("gyration_tensor", tclcommand_analyze_parse_gyration_tensor) + REGISTER_ANALYSIS_WARN("find_principal_axis", tclcommand_analyze_parse_find_principal_axis) REGISTER_ANALYSIS("nbhood", tclcommand_analyze_parse_nbhood); REGISTER_ANALYSIS("distto", tclcommand_analyze_parse_distto); - REGISTER_ANALYSIS("cell_gpb", tclcommand_analyze_parse_cell_gpb); + REGISTER_ANALYSIS_WARN("cell_gpb", tclcommand_analyze_parse_cell_gpb) REGISTER_ANALYSIS("Vkappa", tclcommand_analyze_parse_Vkappa); REGISTER_ANALYSIS("energy", tclcommand_analyze_parse_and_print_energy); REGISTER_ANALYSIS("energy_kinetic", tclcommand_analyze_parse_and_print_energy_kinetic); @@ -2778,18 +2787,18 @@ int tclcommand_analyze(ClientData data, Tcl_Interp *interp, int argc, char **arg // The following analysis commands apply only to the "center of mass" // implementation of virtual sites #ifdef VIRTUAL_SITES_COM - REGISTER_ANALYSIS("energy_kinetic_mol", tclcommand_analyze_parse_and_print_energy_kinetic_mol); - REGISTER_ANALYSIS("pressure_mol", tclcommand_analyze_parse_and_print_pressure_mol); - REGISTER_ANALYSIS("check_mol", tclcommand_analyze_parse_and_print_check_mol); - REGISTER_ANALYSIS("dipmom_mol", tclcommand_analyze_parse_and_print_dipmom_mol); + REGISTER_ANALYSIS_WARN("energy_kinetic_mol", tclcommand_analyze_parse_and_print_energy_kinetic_mol) + REGISTER_ANALYSIS_WARN("pressure_mol", tclcommand_analyze_parse_and_print_pressure_mol) + REGISTER_ANALYSIS_WARN("check_mol", tclcommand_analyze_parse_and_print_check_mol) + REGISTER_ANALYSIS_WARN("dipmom_mol", tclcommand_analyze_parse_and_print_dipmom_mol) #endif #endif REGISTER_ANALYSIS_W_ARG("stress_tensor", tclcommand_analyze_parse_and_print_stress_tensor, 0); REGISTER_ANALYSIS("local_stress_tensor", tclcommand_analyze_parse_local_stress_tensor); - REGISTER_ANALYSIS_W_ARG("p_inst", tclcommand_analyze_parse_and_print_pressure, 1); - REGISTER_ANALYSIS("momentum", tclcommand_analyze_parse_and_print_momentum); - REGISTER_ANALYSIS("bins", tclcommand_analyze_parse_bins); - REGISTER_ANALYSIS("p_IK1", tclcommand_analyze_parse_and_print_p_IK1); + REGISTER_ANALYSIS_W_ARG_WARN("p_inst", tclcommand_analyze_parse_and_print_pressure, 1) + REGISTER_ANALYSIS_WARN("momentum", tclcommand_analyze_parse_and_print_momentum) + REGISTER_ANALYSIS_WARN("bins", tclcommand_analyze_parse_bins) + REGISTER_ANALYSIS_WARN("p_IK1", tclcommand_analyze_parse_and_print_p_IK1) REGISTER_ANALYSIS_W_ARG("re", tclcommand_analyze_parse_re, 0); REGISTER_ANALYSIS_W_ARG("", tclcommand_analyze_parse_re, 1); REGISTER_ANALYSIS_W_ARG("rg", tclcommand_analyze_parse_rg, 0); @@ -2802,27 +2811,27 @@ int tclcommand_analyze(ClientData data, Tcl_Interp *interp, int argc, char **arg REGISTER_ANALYSIS_W_ARG("", tclcommand_analyze_parse_bond_l, 1); REGISTER_ANALYSIS_W_ARG("bond_dist", tclcommand_analyze_parse_bond_dist, 0); REGISTER_ANALYSIS_W_ARG("", tclcommand_analyze_parse_bond_dist, 1); - REGISTER_ANALYSIS_W_ARG("g123", tclcommand_analyze_parse_g123, 1); - REGISTER_ANALYSIS_W_ARG("", tclcommand_analyze_parse_g_av, 1); - REGISTER_ANALYSIS_W_ARG("", tclcommand_analyze_parse_g_av, 2); - REGISTER_ANALYSIS_W_ARG("", tclcommand_analyze_parse_g_av, 3); + REGISTER_ANALYSIS_W_ARG_WARN("g123", tclcommand_analyze_parse_g123, 1) + REGISTER_ANALYSIS_W_ARG_WARN("", tclcommand_analyze_parse_g_av, 1) + REGISTER_ANALYSIS_W_ARG_WARN("", tclcommand_analyze_parse_g_av, 2) + REGISTER_ANALYSIS_W_ARG_WARN("", tclcommand_analyze_parse_g_av, 3) REGISTER_ANALYSIS_W_ARG("formfactor", tclcommand_analyze_parse_formfactor, 0); REGISTER_ANALYSIS_W_ARG("", tclcommand_analyze_parse_formfactor, 1); - REGISTER_ANALYSIS("necklace", tclcommand_analyze_parse_necklace); - REGISTER_ANALYSIS("holes", tclcommand_analyze_parse_holes); - REGISTER_ANALYSIS("distribution", tclcommand_analyze_parse_distribution); - REGISTER_ANALYSIS("vel_distr", tclcommand_analyze_parse_vel_distr); + REGISTER_ANALYSIS_WARN("necklace", tclcommand_analyze_parse_necklace) + REGISTER_ANALYSIS_WARN("holes", tclcommand_analyze_parse_holes) + REGISTER_ANALYSIS_WARN("distribution", tclcommand_analyze_parse_distribution) + REGISTER_ANALYSIS_WARN("vel_distr", tclcommand_analyze_parse_vel_distr) REGISTER_ANALYSIS_W_ARG("rdf", tclcommand_analyze_parse_rdf, 0); REGISTER_ANALYSIS_W_ARG("", tclcommand_analyze_parse_rdf, 1); REGISTER_ANALYSIS_W_ARG("", tclcommand_analyze_parse_rdf, 2); - REGISTER_ANALYSIS("rdfchain", tclcommand_analyze_parse_rdfchain); + REGISTER_ANALYSIS_WARN("rdfchain", tclcommand_analyze_parse_rdfchain) #ifdef ELECTROSTATICS - REGISTER_ANALYSIS("cwvac", tclcommand_analyze_parse_cwvac); + REGISTER_ANALYSIS_WARN("cwvac", tclcommand_analyze_parse_cwvac) #endif - REGISTER_ANALYSIS("structurefactor", tclcommand_analyze_parse_structurefactor); + REGISTER_ANALYSIS_WARN("structurefactor", tclcommand_analyze_parse_structurefactor) REGISTER_ANALYSIS("", tclcommand_analyze_parse_density_profile_av); REGISTER_ANALYSIS("", tclcommand_analyze_parse_diffusion_profile); - REGISTER_ANALYSIS("vanhove", tclcommand_analyze_parse_vanhove); + REGISTER_ANALYSIS_WARN("vanhove", tclcommand_analyze_parse_vanhove) REGISTER_ANALYZE_STORAGE("append", tclcommand_analyze_parse_append); REGISTER_ANALYZE_STORAGE("push", tclcommand_analyze_parse_push); REGISTER_ANALYZE_STORAGE("replace", tclcommand_analyze_parse_replace); @@ -2831,7 +2840,7 @@ int tclcommand_analyze(ClientData data, Tcl_Interp *interp, int argc, char **arg REGISTER_ANALYZE_STORAGE("stored", tclcommand_analyze_parse_stored); REGISTER_ANALYZE_STORAGE("configs", tclcommand_analyze_parse_configs); #ifdef CONFIGTEMP - REGISTER_ANALYSIS("configtemp", tclcommand_analyze_parse_and_print_configtemp); + REGISTER_ANALYSIS_WARN("configtemp", tclcommand_analyze_parse_and_print_configtemp) #endif else { /* the default */ diff --git a/testsuite/Makefile.am b/testsuite/Makefile.am index 70c19b9466a..1203650d74f 100755 --- a/testsuite/Makefile.am +++ b/testsuite/Makefile.am @@ -21,7 +21,9 @@ tests = \ analysis.tcl \ angle.tcl \ bonded_coulomb.tcl \ + collision-detection-angular.tcl \ collision-detection-centers.tcl \ + collision-detection-glue.tcl \ collision-detection-poc.tcl \ comforce.tcl \ comfixed.tcl \ @@ -32,6 +34,7 @@ tests = \ correlation_checkpoint.tcl \ constraints_rhomboid.tcl \ coulomb_cloud_wall.tcl \ + dawaanr-and-dds-gpu.tcl \ dh.tcl \ dielectric.tcl \ dihedral.tcl \ diff --git a/testsuite/collision-detection-angular.tcl b/testsuite/collision-detection-angular.tcl new file mode 100644 index 00000000000..60d47a8d45a --- /dev/null +++ b/testsuite/collision-detection-angular.tcl @@ -0,0 +1,90 @@ +# Copyright (C) 2011,2012,2013,2014 The ESPResSo project +# +# This file is part of ESPResSo. +# +# ESPResSo is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# ESPResSo is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + +# +############################################################# +# # +# Test collision detection with binding of centers of colliding particles +# # +############################################################# +source "tests_common.tcl" + +require_feature "COLLISION_DETECTION" +require_feature "VIRTUAL_SITES_RELATIVE" +require_feature "BOND_ANGLE" + +puts "---------------------------------------------------------------" +puts "- Testcase collision-detection-angular.tcl running on [setmd n_nodes] nodes" +puts "---------------------------------------------------------------" + +# Setup +setmd box_l 10 10 10 + +thermostat off +setmd time_step 0.01 +setmd skin 0 + +# pair bond +inter 2 harmonic 1 1 +# angular bonds +for {set a 0} {$a <= 180} {incr a} { + inter [expr 3 + $a] angle_harmonic 200.0 [expr $a * [PI] / 180] +} + +# triangle around x=0 +part 0 pos -0.366 2 0 +# place close to boundary to check pbc and processor boundaries +part 1 pos 0.5 2.5 0 +part 2 pos 0.5 1.5 0 + +# reverse triangle +part 3 pos 0.366 2 5 +# place close to boundary to check pbc and processor boundaries +part 4 pos -0.5 2.5 5 +part 5 pos -0.5 1.5 5 + +# line +part 6 pos 7 7 7 +part 7 pos 8 7 7 +part 8 pos 9 7 7 + +# Check setting of parameters +setmd min_global_cut 1.0 +on_collision bind_three_particles 1.0 2 3 180 + +set res [on_collision] +if { ! ( ([lindex $res 0] == "bind_three_particles") && (abs([lindex $res 1]-1) <1E-5) + && ([lindex $res 2] == 2) && ([lindex $res 3] == 3) && ([lindex $res 4] == 180)) } { + error_exit "Setting collision_detection parameters for bind_centers does not work" +} + +# Check the actual collision detection +integrate 0 + +for {set p 0} {$p < 9} {incr p} { + puts $p-->[part $p pr bonds] +} + +puts "-------------------------------------------------------------------" +# Integrate again and make sure, no extra bonds are added +integrate 0 recalc_forces + +for {set p 0} {$p < 9} {incr p} { + puts $p-->[part $p pr bonds] +} + +exit 0 diff --git a/testsuite/collision-detection-centers.tcl b/testsuite/collision-detection-centers.tcl index 90332d43f6f..22f64370790 100755 --- a/testsuite/collision-detection-centers.tcl +++ b/testsuite/collision-detection-centers.tcl @@ -1,24 +1,24 @@ # Copyright (C) 2011,2012,2013,2014 The ESPResSo project # -# This file is part of ESPResSo. +# this file is part of espresso. # -# ESPResSo is free software: you can redistribute it and/or modify -# it under the terms of the GNU General Public License as published by -# the Free Software Foundation, either version 3 of the License, or +# espresso is free software: you can redistribute it and/or modify +# it under the terms of the gnu general public license as published by +# the free software foundation, either version 3 of the license, or # (at your option) any later version. # -# ESPResSo is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. +# espresso is distributed in the hope that it will be useful, +# but without any warranty; without even the implied warranty of +# merchantability or fitness for a particular purpose. see the +# gnu general public license for more details. # -# You should have received a copy of the GNU General Public License -# along with this program. If not, see . +# you should have received a copy of the gnu general public license +# along with this program. if not, see . # ############################################################# # # -# Test collision detection with binding of centers of colliding particles +# test collision detection with binding of centers of colliding particles # # ############################################################# source "tests_common.tcl" @@ -32,14 +32,14 @@ puts "---------------------------------------------------------------" puts "- Testcase collision-detection-centers.tcl running on [setmd n_nodes] nodes" puts "---------------------------------------------------------------" -# Setup +# setup setmd box_l 10 10 10 thermostat off setmd time_step 0.01 inter 3 harmonic 2 2 inter 0 0 lennard-jones 0.0001 2 2.1 auto -inter 7 harmonic 1 1 +inter 7 harmonic 2 1 setmd skin 0 part 0 pos 9 0 0 # place close to boundary to check pbc and processor boundaries @@ -68,35 +68,35 @@ proc analyze_topology {bond_type {check_others 0}} { return [lsort $bonded] } -# Test default setting +# test default setting if { "[on_collision]" != "off" } { - error_exit "Collision detection should be off by default." + error_exit "collision detection should be off by default." } -# Test switching it off +# test switching it off on_collision off if { "[on_collision]" != "off" } { - error_exit "Disabling collision_detection does not work" + error_exit "disabling collision_detection does not work" } -# Make sure, it doesn't do anything when turned off +# make sure, it doesn't do anything when turned off integrate 0 set bonds [analyze_topology "" 1] if {$bonds != ""} { - error_exit "Bonds were created when collision detection was off." + error_exit "bonds were created when collision detection was off." } -# Check setting of parameters +# check setting of parameters setmd min_global_cut 1.0 -on_collision bind_centers 2.0 7 +on_collision bind_centers 2.0 7 set res [on_collision] -if { ! ( ([lindex $res 0] == "bind_centers") && (abs([lindex $res 1]-2) <1E-5) && ([lindex $res 2] == 7)) } { - error_exit "Setting collision_detection parameters for bind_centers does not work" +if { ! ( ([lindex $res 0] == "bind_centers") && (abs([lindex $res 1]-2) <1e-5) && ([lindex $res 2] == 7)) } { + error_exit "setting collision_detection parameters for bind_centers does not work" } -# Check the actual collision detection +# check the actual collision detection integrate 0 # Check, whether the bonds are correct @@ -120,7 +120,7 @@ part 0 pos 0.5 0 0 integrate 0 -# Check, whether the bonds are still correct, not doubled +# check, whether the bonds are still correct, not doubled set bonds [analyze_topology 7 1] if {$bonds != "{0 1}"} { error_exit "bond double on exchange: bonds are $bonds" @@ -137,7 +137,7 @@ if {![catch {integrate 0} err]} { set bonds "" foreach exception [lrange $err 1 end] { - if {[regexp {collision between particles (\d+) and (\d+)} $exception -> id1 id2]} { + if {[regexp { ERROR: collision between particles (\d+) and (\d+)} $exception -> id1 id2]} { lappend bonds "$id1 $id2" } else { error_exit "unexpected exception $exception" diff --git a/testsuite/collision-detection-glue.tcl b/testsuite/collision-detection-glue.tcl new file mode 100644 index 00000000000..3fc3c772126 --- /dev/null +++ b/testsuite/collision-detection-glue.tcl @@ -0,0 +1,186 @@ +# Copyright (C) 2011,2012 The ESPResSo project +# +# This file is part of ESPResSo. +# +# ESPResSo is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# ESPResSo is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + +# +############################################################# +# # +# Test collision detection in the "glue to surface>" mode +# # +############################################################# +source "tests_common.tcl" + +require_feature "VIRTUAL_SITES_RELATIVE" +require_feature "COLLISION_DETECTION" +require_feature "ADRESS" off +require_max_nodes_per_side {1 1 1} + +puts "---------------------------------------------------------------" +puts "- Testcase collision-detection-glue.tcl running on 1 nodes" +puts "---------------------------------------------------------------" + +# Setup +setmd box_l 19 19 19 + +thermostat off +setmd time_step 0.01 +inter 2 harmonic 1 1 +inter 3 harmonic 1 0.0001 +setmd skin 0 +part 0 pos 0 0 0 type 2 +part 1 pos 5.5 0 0 type 3 +part 2 pos 3 0 0 type 0 +inter 2 3 lennard-jones 0.001 5.6 5.7 auto +# analyze the bonding structure for pair bonds +proc analyze_topology {bond_type {check_others 0}} { + set bonded "" + for {set i 0} {$i <= [setmd max_part]} {incr i} { + foreach bond [lindex [part $i pr bond] 0] { + if {[lindex $bond 0] == $bond_type} { + set j [lindex $bond 1] + if {$i < $j} { + lappend bonded "$i $j" + } { + lappend bonded "$j $i" + } + } { + if {$check_others} { + error_exit "bond $bond at particle $i of unexpected type found" + } + } + } + } + return [lsort $bonded] +} + +# Test default setting +if { "[on_collision]" != "off" } { + error_exit "Collision detection should be off by default." +} + +# Test switching it off +on_collision off +if { "[on_collision]" != "off" } { + error_exit "Disabling collision_detection does not work" +} + +# Make sure, it doesn't do anything when turned off +integrate 0 + +set bonds [analyze_topology "" 1] +if {$bonds != ""} { + error_exit "Bonds were created when collision detection was off." +} + +# Check setting of parameters +setmd min_global_cut 5.0 +on_collision glue_to_surface 5.5 2 3 1 2 3 4 1.0 + +set res [on_collision] + +if { ! ($res == "glue_to_surface 5.500000 2 3 1 2 3 4 1.000000") } { + error_exit "Setting collision_detection parameters does not work: $res" +} + +# Check the actual collision detection +integrate 0 recalc_forces + +# Check bonds between colliding particles +set bonds [analyze_topology 2] +if {$bonds != "{0 1}"} { + error_exit "bond between colliding particles not created as it should: bonds are $bonds" +} + +# Check bond between the glued particle and the virtual site +set bonds [analyze_topology 3] +if {$bonds != "{0 3}"} { + error_exit "bond between glued particle and vs not created as it should: bonds are $bonds" +} + +# Check if the virtual site has the correct settings +if { [part 3 print virtual] != 1 } { + error_exit "The supposed virtual particle doesn't have the virtual flag set." +} +# Is the vs attached correctly +set vs_info [lrange [part 3 print vs_relative] 0 1] +if {$vs_info != "1 4.500000" } { + error_exit "the vs_relative params are wrong: $vs_info" +} + +# Integrate again and make sure, no extra bonds are added +# enforce force recalculation +integrate 0 recalc_forces + +# Check, whether the bonds are still correct, not doubled +set bonds [analyze_topology 2] +if {$bonds != "{0 1}"} { + error_exit "After 2nd run: bond between colliding particles not created as it should: bonds are $bonds" +} + +# Check bond between the glued particle and the virtual site +set bonds [analyze_topology 3] +if {$bonds != "{0 3}"} { + error_exit "After 2nd run: bond between glued particle and vs not created as it should: bonds are $bonds" +} + + + +# Check whether the number of particles is correct (3 normal +1 vs =4) +if {[setmd n_part] != 4} { + error_exit "Incorrect number of particles [setmd n_part] in the simulation. Too many or too few virtual sites were created." +} + +# Check the particle type of the virtual sites +if { ! ([part 3 print type]==1) } { + error_exit "type of virtual sites is incorrect." +} +if { ! ([part 0 print type]==4) } { + error_exit "type of glued particle is incorrect: [part 0 print type]." +} + +# test exception, generating another collision +part 2 pos 11 0 0 type 2 +on_collision exception glue_to_surface 5.5 2 3 1 2 3 4 1.0 +if {![catch {integrate 1} err]} { + error_exit "no exception was thrown at collision, although requested" +} + +set bonds "" +foreach exception [lrange $err 1 end] { + if {[lrange $exception 1 3] != "collision between particles"} { + error_exit "unexpected exception $exception" + } + lappend bonds "[lindex $exception 4] [lindex $exception 6]" +} +set bonds [lsort $bonds] + +if {$bonds != "{1 2}"} { + error_exit "exception bonds $bonds wrong" +} + +# Check, whether the bonds are also correct +# Between centers of colliding particles +set bonds [analyze_topology 2] +if {$bonds != "{0 1} {1 2}"} { + error_exit "bonds not correctly created: bonds are $bonds" +} +# Between glued particle and vs +set bonds [analyze_topology 3] +if {$bonds != "{0 3} {2 4}"} { + error_exit "bonds not correctly created: bonds are $bonds" +} + +exit 0 diff --git a/testsuite/collision-detection-poc.tcl b/testsuite/collision-detection-poc.tcl index 1e7619c1aa2..cde8e655cc2 100755 --- a/testsuite/collision-detection-poc.tcl +++ b/testsuite/collision-detection-poc.tcl @@ -88,13 +88,14 @@ if {$bonds != ""} { setmd min_global_cut 1.0 on_collision bind_at_point_of_collision 1.0 2 3 1 + set res [on_collision] if { ! ( ([lindex $res 0] == "bind_at_point_of_collision") && (abs([lindex $res 1]-1) <1E-5) && ([lindex $res 2] == 2) && ([lindex $res 3] == 3) && ([lindex $res 4] == 1)) } { error_exit "Setting collision_detection parameters for bind_centers does not work" } # Check the actual collision detection -integrate 0 +integrate 0 recalc_forces # Check, whether the bonds are correct. No strict checking, there are also the # virtual particles @@ -130,6 +131,37 @@ if { (! ([part 3 print type]==1 && [part 4 print type]==1))} { error_exit "type of virtual sites is incorrect." } + +# Check position and vs_relative settings of virtual sites +set n_related_to_0 0 +set n_related_to_1 0 +for {set i 3} {$i <=4} {incr i} { + set vs_r [part $i print vs_relative] + set relto [lindex $vs_r 0] + set dist [lindex $vs_r 1] + + if { abs($dist -0.5)>1E-5 } { + error_exit "Distance between vs particle $i and particle $relto wrong: $dist" + } + if { $relto == 1 } { + incr n_related_to_1 + } else { + if { $relto == 0 } { + incr n_related_to_0 + } else { + error_exit "Vs $i should not be relatex to $relto" + } + } +} +if { ($n_related_to_0 != 1) || ($n_related_to_1 != 1) } { + error_exit "Exactly one vs is supposed to be related to part 0 and one to part 1." +} + + + + + + # test exception, generating another collision part 2 pos 2 0 0 on_collision exception bind_at_point_of_collision 1.0 2 3 1 diff --git a/testsuite/dawaanr-and-dds-gpu.tcl b/testsuite/dawaanr-and-dds-gpu.tcl new file mode 100644 index 00000000000..cf021bdcb53 --- /dev/null +++ b/testsuite/dawaanr-and-dds-gpu.tcl @@ -0,0 +1,168 @@ + + +# Copyright (C) 2010,2011,2012,2013,2014 The ESPResSo project +# Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010 +# Max-Planck-Institute for Polymer Research, Theory Group +# +# This file is part of ESPResSo. +# +# ESPResSo is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# ESPResSo is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# + +proc stopAll {} { + for {set i 0} {$i <=[setmd max_part] } {incr i} { + if { [part $i print pos] != "na" } { + part $i v 0 0 0 omega 0 0 0 + } + } +} + +proc error { msg } { + puts $msg +} + +source "tests_common.tcl" +require_feature "DIPOLES" +require_feature "CUDA" +require_feature "PARTIAL_PERIODIC" +require_max_nodes_per_side 1 + +source "tests_common.tcl" + +set tcl_precision 14 + +proc vectorsTheSame {a b} { + set tol 3E-3 + set diff [vecsub $a $b] + if { [veclen $diff] > $tol } { + return 0 + puts "Difference: [veclen $diff]" + } + return 1 +} + + +puts "----------------------------------------------" +puts "- Testcase dawaanr-and-dds-gpu.tcl" +puts "----------------------------------------------" + +set pf_dds_gpu 2.34 +set pf_dawaanr 3.524 +set ratio_dawaanr_dds_gpu [expr $pf_dawaanr / $pf_dds_gpu] + +set l 15 +setmd box_l $l $l $l +setmd periodic 0 0 0 +setmd time_step 0.0001 +setmd skin 0.1 + + +foreach n { 110 111 540 541 } { + + + +puts "$n particles " +set dipole_modulus 1.3 +for {set i 0 } {$i < $n} {incr i} { + set posx [expr $l*[expr [t_random]]] + set posy [expr $l*[expr [t_random]]] + set posz [expr $l*[expr [t_random]]] + set costheta [expr 2*[expr [t_random]] - 1] + set sintheta [expr sin(acos($costheta))] + set phi [expr 2*3.1415926536*[expr [t_random]]] + set dipx [expr $sintheta*cos($phi)*$dipole_modulus] + set dipy [expr $sintheta*sin($phi)*$dipole_modulus] + set dipz [expr $costheta*$dipole_modulus] + + part [expr 2*$i] pos $posx $posy $posz type 0 dip $dipx $dipy $dipz v 0 0 0 +} +inter 0 0 lennard-jones 10 0.5 0.55 auto +thermostat langevin 0 10 +minimize_energy 0 500 0.1 0.1 +stopAll + +inter 0 0 lennard-jones 0 0 0 0 + + + +setmd skin 0 +setmd time_step 0.01 + +inter magnetic $pf_dawaanr dawaanr +integrate 0 recalc_forces + +set dawaanr_f "" +set dawaanr_t "" + +for {set i 0} {$i<$n} {incr i} { + lappend dawaanr_f [part [expr 2*$i] print f] + lappend dawaanr_t [part [expr 2*$i] print torque] +} +set dawaanr_e [analyze energy total] + +inter magnetic 0 + +integrate 0 recalc_forces + +inter magnetic $pf_dds_gpu dds-gpu + +integrate 0 recalc_forces + +set ddsgpu_f "" +set ddsgpu_t "" + +for {set i 0} {$i<$n} {incr i} { + lappend ddsgpu_f [part [expr 2*$i] print f] + lappend ddsgpu_t [part [expr 2*$i] print torque] +} + +set ddsgpu_e [analyze energy total] + +# compare + +for {set i 0} {$i<$n} {incr i} { + if { ! [vectorsTheSame [lindex $dawaanr_f $i] [vecscale $ratio_dawaanr_dds_gpu [lindex $ddsgpu_f $i]]] } { + error "Forces on particle don't match. $i [vecsub [lindex $dawaanr_f $i] [vecscale $ratio_dawaanr_dds_gpu [lindex $ddsgpu_f $i]]]\n[lindex $dawaanr_f $i] [vecscale $ratio_dawaanr_dds_gpu [lindex $ddsgpu_f $i]]" + } + if { ! [vectorsTheSame [lindex $dawaanr_t $i] [vecscale $ratio_dawaanr_dds_gpu [lindex $ddsgpu_t $i]]] } { + error "torques on particle $i don't match. [vecsub [lindex $dawaanr_t $i] [lindex $ddsgpu_t $i]] " + } +} + +if { abs($dawaanr_e - $ddsgpu_e*$ratio_dawaanr_dds_gpu) > 0.001 } { + error "Energies for dawaanr $dawaanr_e and dds_gpu $ddsgpu_e don't match." +} + +# Does the energy stay constant when swapping particles +for {set i 0} {$i<1000} {incr i} { + set a [expr int(rand() *$n)*2] + set b [expr int(rand() *$n)*2] + set dipa [part $a print dip] + set posa [part $a print pos] + eval "part $a pos [part $b print pos] dip [part $b print dip]" + eval "part $b pos $posa dip $dipa" + integrate 0 recalc_forces + set E [analyze energy total] + if { abs($E -$ddsgpu_e) > 1E-3 } { + error "energy mismatch: $E $ddsgpu_e" + } +} + + +integrate 0 recalc_forces +inter magnetic 0 +part deleteall + + +}