From f04fb8edafcfe6c960f9a28cf9b934812bbf1561 Mon Sep 17 00:00:00 2001 From: Rudolf Weeber Date: Wed, 29 Dec 2010 15:53:41 +0100 Subject: [PATCH 1/8] Make rotation switchable on a per-particle level When SWITCHABLE_ROTATION is defined in myconfig.h, part ... rotation 1 activates rotation for a specific particle. SWITCHABLE_ROTATION is incompatible with pressure calculation, as there is no method to count the number of degrees of freedom, yet. --- communication.c | 41 +++++++++++++++++++++++-- communication.h | 13 +++++++- doc/ug/features.tex | 1 + doc/ug/part.tex | 3 ++ energy.h | 5 +++ myconfig-sample.h | 2 ++ particle_data.c | 74 +++++++++++++++++++++++++++++++++++++++++++++ particle_data.h | 6 ++++ pressure.c | 5 +++ pressure.h | 3 ++ rotation.c | 12 ++++++++ 11 files changed, 162 insertions(+), 3 deletions(-) diff --git a/communication.c b/communication.c index 4278ce16a52..86ad3866b8d 100644 --- a/communication.c +++ b/communication.c @@ -191,9 +191,11 @@ typedef void (SlaveCallback)(int node, int param); #define REQ_SET_RINERTIA 54 /** Action number for \ref mpi_send_mu_E. */ #define REQ_SET_MU_E 55 +/** Action number for \ref mpi_send_rotation. */ +#define REQ_SET_ROTATION 56 /** Total number of action numbers. */ -#define REQ_MAXIMUM 56 +#define REQ_MAXIMUM 57 /*@}*/ @@ -318,7 +320,8 @@ static SlaveCallback *slave_callbacks[] = { mpi_iccp3m_iteration_slave, /* 52: REQ_ICCP3M_ITERATION */ mpi_iccp3m_init_slave, /* 53: REQ_ICCP3M_INIT */ mpi_send_rotational_inertia_slave,/* 54: REQ_SET_RINERTIA */ - mpi_send_mu_E_slave, /* 55: REQ_SET_MU_E */ + mpi_send_mu_E_slave, /* 55: REQ_SET_MU_E */ + mpi_send_rotation, /* 56: REQ_SET_ROTTION */ }; /** Names to be printed when communication debugging is on. */ @@ -389,6 +392,7 @@ char *names[] = { "REQ_ICCP3M_INIT", /* 53 */ "SET_RINERTIA", /* 54 */ "SET_MU_E", /* 55 */ + "SET_ROTATION", /* 56 */ }; /** the requests are compiled here. So after a crash you get the last issued request */ @@ -1125,6 +1129,39 @@ void mpi_send_virtual_slave(int pnode, int part) #endif } +// ******************************** + +void mpi_send_rotation(int pnode, int part, int rot) +{ +#ifdef SWITCHABLE_ROTATION + mpi_issue(REQ_SET_ROTATION, pnode, part); + + if (pnode == this_node) { + Particle *p = local_particles[part]; + p->p.rotation = rot; + } + else { + MPI_Send(&rot, 1, MPI_INT, pnode, REQ_SET_ROTATION, MPI_COMM_WORLD); + } + + on_particle_change(); +#endif +} + +void mpi_send_rotation_slave(int pnode, int part) +{ +#ifdef SWITCHABLE_ROTATION + if (pnode == this_node) { + Particle *p = local_particles[part]; + MPI_Status status; + MPI_Recv(&p->p.rotation, 1, MPI_INT, 0, REQ_SET_ROTATION, + MPI_COMM_WORLD, &status); + } + + on_particle_change(); +#endif +} + /********************* REQ_SET_BOND ********/ int mpi_send_bond(int pnode, int part, int *bond, int delete) { diff --git a/communication.h b/communication.h index 9a1a168cf4f..fcb61938e4d 100644 --- a/communication.h +++ b/communication.h @@ -180,7 +180,17 @@ void mpi_send_rotational_inertia(int node, int part, double rinertia[3]); */ void mpi_send_quat(int node, int part, double quat[4]); -/** Issue REQ_SET_LAMBDA: send particle angular velocity. + +/** Issue REQ_SET_ROTATION: send particle rotation flag + Also calls \ref on_particle_change. + \param part the particle. + \param node the node it is attached to. + \param rot the rotation flag +*/ +void mpi_send_rotation(int pnode, int part, int rot); + + +/* Issue REQ_SET_LAMBDA: send particle angular velocity. Also calls \ref on_particle_change. \param part the particle. \param node the node it is attached to. @@ -197,6 +207,7 @@ void mpi_send_omega(int node, int part, double omega[3]); void mpi_send_torque(int node, int part, double torque[3]); #endif + #ifdef DIPOLES /** Issue REQ_SET_DIP: send particle dipole orientation. Also calls \ref on_particle_change. diff --git a/doc/ug/features.tex b/doc/ug/features.tex index 60bacc6987e..f1422b11edc 100644 --- a/doc/ug/features.tex +++ b/doc/ug/features.tex @@ -60,6 +60,7 @@ \section{General features} freedom, which for example means that the kinetic energy changes at constant temperature is twice as large. \todo{Docs for rotation missing} +\item \newfeature{SWITCHABLE_ROTATION} Makes the integration of rotational degrees of freedom switchable on a per-particle level. \item \newfeature{EXTERNAL\_FORCES} Allows to define an arbitrary constant force for each particle individually. Also allows to fix individual coordinates of particles, \eg keep them at a fixed diff --git a/doc/ug/part.tex b/doc/ug/part.tex index 19470b6ff8e..9ef28b14ad9 100644 --- a/doc/ug/part.tex +++ b/doc/ug/part.tex @@ -44,6 +44,7 @@ \subsection{Defining particle properties} \require{5}{\opt{mass \var{mass}}} \require{6}{\opt{dipm \var{moment}}} \require{6}{\opt{dip \var{dx} \var{dy} \var{dz}}} + \require{7}{\opt{rotation \var{rot}}} \begin{features} \required[1]{ELECTROSTATICS} \required[2]{ROTATION} @@ -51,6 +52,7 @@ \subsection{Defining particle properties} \required[4]{EXCLUSION} \required[5]{MASS} \required[6]{DIPOLES} + \required[6]{SWITCHABLE_ROTATION} \end{features} \end{essyntax} @@ -113,6 +115,7 @@ \subsection{Defining particle properties} \item[\opt{dipm \var{moment}}] Sets the dipol moment of this particle to $moment$. \item[\opt{dip \var{dx} \var{dy} \var{dz}}] Sets the orientation of the dipol axis to $(dx,dy,dz)$. + \item \opt{rotation \var{rot}} specifies whether a particle's rotational degrees of freedom are integrated (value of 1) or not (0). If set to zero, the content of the torque and omega variables are meaningless. \end{arguments} \subsection{Getting particle properties} diff --git a/energy.h b/energy.h index fc890934bce..1bca54b8424 100644 --- a/energy.h +++ b/energy.h @@ -418,6 +418,10 @@ MDINLINE void add_kinetic_energy(Particle *p1) energy.data.e[0] += (SQR(p1->m.v[0]) + SQR(p1->m.v[1]) + SQR(p1->m.v[2]))*PMASS(*p1); #ifdef ROTATION +#ifdef SWITCHABLE_ROTATION +if (p1->p.rotation) +#endif +{ #ifdef ROTATIONAL_INERTIA /* the rotational part is added to the total kinetic energy; Here we use the rotational inertia */ @@ -430,6 +434,7 @@ MDINLINE void add_kinetic_energy(Particle *p1) at the moment, we assume unit inertia tensor I=(1,1,1) */ energy.data.e[0] += (SQR(p1->m.omega[0]) + SQR(p1->m.omega[1]) + SQR(p1->m.omega[2]))*time_step*time_step; #endif +} #endif } diff --git a/myconfig-sample.h b/myconfig-sample.h index 28443dab44b..01ec317eddd 100644 --- a/myconfig-sample.h +++ b/myconfig-sample.h @@ -24,6 +24,8 @@ /* #define PARTIAL_PERIODIC */ /* #define ELECTROSTATICS */ /* #define ROTATION */ +// Enable the following, if you want to enable rotation on a per-particle level +// $define SWITCHABLE_ROTATION /* #define ROTATIONAL_INTERIA */ /* #define DIPOLES */ /* #define MAGNETOSTATICS */ diff --git a/particle_data.c b/particle_data.c index 52a3417b9f2..7d03b28bca6 100644 --- a/particle_data.c +++ b/particle_data.c @@ -109,6 +109,10 @@ void init_particle(Particle *part) part->p.rinertia[1] = 1.0; part->p.rinertia[2] = 1.0; #endif +#ifdef SWITCHABLE_ROTATION + part->p.rotation =0; +#endif + #ifdef ELECTROSTATICS part->p.q = 0.0; @@ -545,6 +549,14 @@ void tclcommand_part_print_virtual(Particle *part, char *buffer, Tcl_Interp *int } #endif +#ifdef SWITCHABLE_ROTATION +void tclcommand_part_print_rotation(Particle *part, char *buffer, Tcl_Interp *interp) +{ + sprintf(buffer,"%i", part->p.rotation); + Tcl_AppendResult(interp, buffer, (char *)NULL); +} +#endif + void tclcommand_part_print_v(Particle *part, char *buffer, Tcl_Interp *interp) { /* unscale velocities ! */ @@ -823,6 +835,13 @@ int tclprint_to_result_Particle(Tcl_Interp *interp, int part_num) tclcommand_part_print_torque(&part, buffer, interp); #endif +#ifdef SWITCHABLE_ROTATION + Tcl_AppendResult(interp, " rotation ", (char *)NULL); + tclcommand_part_print_rotation(&part, buffer, interp); +#endif + + + #ifdef ROTATIONAL_INERTIA /* print information about rotational inertia */ Tcl_AppendResult(interp, " rinertia ", (char *)NULL); @@ -971,6 +990,10 @@ int tclcommand_part_parse_print(Tcl_Interp *interp, int argc, char **argv, else if (ARG0_IS_S("tbf")) tclcommand_part_print_torque_body_frame(&part, buffer, interp); #endif +#ifdef SWITCHABLE_ROTATION + else if (ARG0_IS_S("rotation")) + tclcommand_part_print_rotation(&part, buffer, interp); +#endif #ifdef ROTATIONAL_INERTIA else if (ARG0_IS_S("rinertia")) @@ -1136,6 +1159,33 @@ int tclcommand_part_parse_rotational_inertia(Tcl_Interp *interp, int argc, char } #endif +#ifdef SWITCHABLE_ROTATION +int tclcommand_part_parse_rotation(Tcl_Interp *interp, int argc, char **argv, + int part_num, int * change) +{ + int rot; + + *change = 1; + + if (argc < 1) { + Tcl_AppendResult(interp, "rotation requires 1 argument", (char *) NULL); + return TCL_ERROR; + } + + /* set rotation flag */ + if (! ARG0_IS_I(rot)) + return TCL_ERROR; + + if (set_particle_rotation(part_num, rot) == TCL_ERROR) { + Tcl_AppendResult(interp, "set particle position first", (char *)NULL); + + return TCL_ERROR; + } + + return TCL_OK; +} +#endif + #ifdef DIPOLES int tclcommand_part_parse_dipm(Tcl_Interp *interp, int argc, char **argv, int part_num, int * change) @@ -1897,6 +1947,11 @@ int tclcommand_part_parse_cmd(Tcl_Interp *interp, int argc, char **argv, err = tclcommand_part_parse_rotational_inertia(interp, argc-1, argv+1, part_num, &change); #endif +#ifdef SWITCHABLE_ROTATION + else if (ARG0_IS_S("rotation")) + err = tclcommand_part_parse_rotation(interp, argc-1, argv+1, part_num, &change); +#endif + #ifdef DIPOLES else if (ARG0_IS_S("dip")) @@ -2139,6 +2194,25 @@ int set_particle_rotational_inertia(int part, double rinertia[3]) } #endif + +#ifdef SWITCHABLE_ROTATION +int set_particle_rotation(int part, int rot) +{ + int pnode; + if (!particle_node) + build_particle_node(); + + if (part < 0 || part > max_seen_particle) + return TCL_ERROR; + pnode = particle_node[part]; + + if (pnode == -1) + return TCL_ERROR; + mpi_send_rotation(pnode, part, rot); + return TCL_OK; +} +#endif + #ifdef DIPOLES int set_particle_dipm(int part, double dipm) { diff --git a/particle_data.h b/particle_data.h index ff086a0bfe1..7b6b58a4ae6 100644 --- a/particle_data.h +++ b/particle_data.h @@ -78,6 +78,12 @@ typedef struct { double rinertia[3]; #endif +#ifdef SWITCHABLE_ROTATION + // Determines, wether a particle's rotational degrees of freedom are + // integrated + int rotation; +#endif + #ifdef ELECTROSTATICS /** charge. */ double q; diff --git a/pressure.c b/pressure.c index 30d34a092f6..97c7cfb91d1 100644 --- a/pressure.c +++ b/pressure.c @@ -131,6 +131,11 @@ void pressure_calc(double *result, double *result_t, double *result_nb, double * } /* rescale kinetic energy (=ideal contribution) */ #ifdef ROTATION +#ifdef SWITCHABLE_ROTATION + fprintf(stderr, "Switching rotation on the particle level (#define SWITCHABLE_ROTATION) and pressure calculation are incompatible.\n"); +#endif + + virials.data.e[0] /= (6.0*volume*time_step*time_step); #else virials.data.e[0] /= (3.0*volume*time_step*time_step); diff --git a/pressure.h b/pressure.h index e30c198f521..346078de878 100644 --- a/pressure.h +++ b/pressure.h @@ -525,6 +525,9 @@ MDINLINE void add_kinetic_virials(Particle *p1,int v_comp) virials.data.e[0] += (SQR(p1->m.v[0]) + SQR(p1->m.v[1]) + SQR(p1->m.v[2]))*PMASS(*p1); #ifdef ROTATION +#ifdef SWITCHABLE_ROTATION +if (p1->p.rotation) +#endif virials.data.e[0] += (SQR(p1->m.omega[0]) + SQR(p1->m.omega[1]) + SQR(p1->m.omega[2]))*SQR(time_step); #endif diff --git a/rotation.c b/rotation.c index a25b3b1789d..d33a3f56ebc 100644 --- a/rotation.c +++ b/rotation.c @@ -239,6 +239,10 @@ void propagate_omega_quat() p = cell->part; np = cell->n; for(i = 0; i < np; i++) { +#ifdef SWITCHABLE_ROTATION + if (!p[i].p.rotation) + continue; +#endif double Qd[4], Qdd[4], S[3], Wd[3]; define_Qdd(&p[i], Qd, Qdd, S, Wd); @@ -280,6 +284,10 @@ void convert_torqes_propagate_omega() p = cell->part; np = cell->n; for(i = 0; i < np; i++) { +#ifdef SWITCHABLE_ROTATION + if (!p[i].p.rotation) + continue; +#endif double A[9]; define_rotation_matrix(&p[i], A); @@ -351,6 +359,10 @@ void convert_initial_torques() p = cell->part; np = cell->n; for(i = 0; i < np; i++) { +#ifdef SWITCHABLE_ROTATION + if (!p[i].p.rotation) + continue; +#endif double A[9]; define_rotation_matrix(&p[i], A); From 2ec2d57e6e02585c9e85d38a28f58761afa339a4 Mon Sep 17 00:00:00 2001 From: Rudolf Weeber Date: Fri, 31 Dec 2010 09:21:31 +0100 Subject: [PATCH 2/8] SWITCHABLE_ROTATION: Added to feature list in config.c --- config.c | 3 +++ 1 file changed, 3 insertions(+) diff --git a/config.c b/config.c index 288d71b0020..ef6ee0cab2b 100644 --- a/config.c +++ b/config.c @@ -134,6 +134,9 @@ int tclcallback_compilation(Tcl_Interp *interp) #ifdef ROTATION Tcl_AppendResult(interp, "{ ROTATION } ", (char *) NULL); #endif +#ifdef SWITCHABLE_ROTATION + Tcl_AppendResult(interp, "{ SWITCHABLE_ROTATION } ", (char *) NULL); +#endif #ifdef DIPOLES Tcl_AppendResult(interp, "{ DIPOLES } ", (char *) NULL); #endif From dee673f01b7a2bd3f16a61e9782274a301c91ce7 Mon Sep 17 00:00:00 2001 From: Rudolf Weeber Date: Wed, 26 Jan 2011 16:24:27 +0100 Subject: [PATCH 3/8] Dipolar all with all - no replicas: Respect boundary conditions set by setmd --- magnetic_non_p3m__methods.c | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/magnetic_non_p3m__methods.c b/magnetic_non_p3m__methods.c index 95f795c215f..7d7e5061516 100644 --- a/magnetic_non_p3m__methods.c +++ b/magnetic_non_p3m__methods.c @@ -105,13 +105,6 @@ int DAWAANR_sanity_checks() char *errtxt; #endif -#ifdef PARTIAL_PERIODIC - if (!PERIODIC(0) || !PERIODIC(1) || !PERIODIC(2)) { - errtxt = runtime_error(128); - ERROR_SPRINTF(errtxt, "{006 DAWAANR requires periodicity 1 1 1} "); - return 1; - } -#endif return 0; } @@ -202,10 +195,20 @@ double dawaanr_calculations(int force_flag, int energy_flag) { u=0; for(i=0;i Date: Tue, 22 Feb 2011 12:13:02 +0100 Subject: [PATCH 4/8] SWITCHABLE_ROTATION implies ROTATION. Conflicts: communication.c config.h --- communication.c | 3 ++- config.h | 18 ++++++++++++++++++ 2 files changed, 20 insertions(+), 1 deletion(-) diff --git a/communication.c b/communication.c index 86ad3866b8d..f9353852291 100644 --- a/communication.c +++ b/communication.c @@ -260,6 +260,7 @@ void mpi_ljangle_cap_forces_slave(int node, int parm); void mpi_send_virtual_slave(int node, int parm); void mpi_bcast_tf_params_slave(int node, int parm); void mpi_send_rotational_inertia_slave(int node, int parm); +void mpi_send_rotation_slave(int node, int parm); /*@}*/ /** A list of which function has to be called for @@ -321,7 +322,7 @@ static SlaveCallback *slave_callbacks[] = { mpi_iccp3m_init_slave, /* 53: REQ_ICCP3M_INIT */ mpi_send_rotational_inertia_slave,/* 54: REQ_SET_RINERTIA */ mpi_send_mu_E_slave, /* 55: REQ_SET_MU_E */ - mpi_send_rotation, /* 56: REQ_SET_ROTTION */ + mpi_send_rotation_slave, /* 56: REQ_SET_ROTATION */ }; /** Names to be printed when communication debugging is on. */ diff --git a/config.h b/config.h index ac2f70c7dcf..60a28f3cb33 100644 --- a/config.h +++ b/config.h @@ -189,6 +189,24 @@ #define CONSTRAINTS #endif +#if defined(VIRTUAL_SITES_COM) || defined(VIRTUAL_SITES_RELATIVE) +#define VIRTUAL_SITES +#endif + +#ifdef VIRTUAL_SITES_RELATIVE +#ifndef ROTATION +#define ROTATION +#endif +#endif + +#ifdef SWITCHABLE_ROTATION +#ifndef ROTATION +#define ROTATION +#endif +#endif + +/*@}*/ + /********************************************/ /* \name exported functions of config.c */ /********************************************/ From fff407cdbc60eecd229a661523abca6caf4d7550 Mon Sep 17 00:00:00 2001 From: Rudolf Weeber Date: Mon, 9 Jan 2012 16:09:49 +0100 Subject: [PATCH 5/8] in add_constraint_force, exit early if no constraints are active --- src/constraint.h | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/constraint.h b/src/constraint.h index 81c521fdfdf..9d1054a470a 100644 --- a/src/constraint.h +++ b/src/constraint.h @@ -1174,6 +1174,8 @@ MDINLINE void reflect_particle(Particle *p1, double *distance_vec, int reflectin MDINLINE void add_constraints_forces(Particle *p1) { + if (n_constraints==0) + return; int n, j; double dist, vec[3], force[3], torque1[3], torque2[3]; From 5768b9961edfd88bfbbe0f355ecf0b48cbc19bf5 Mon Sep 17 00:00:00 2001 From: Rudolf Weeber Date: Mon, 9 Jan 2012 16:11:51 +0100 Subject: [PATCH 6/8] The code from check_forces() which checks if a force is nan into rescale_forces...() to avoid an additional iteration over all particles --- src/forces.c | 23 ++--------------------- src/forces.h | 21 +++++++++++++++++++++ src/integrate.c | 2 ++ 3 files changed, 25 insertions(+), 21 deletions(-) diff --git a/src/forces.c b/src/forces.c index 0f582a43bcd..3c25d98e1af 100644 --- a/src/forces.c +++ b/src/forces.c @@ -115,7 +115,6 @@ void force_calc() #ifdef COMFIXED calc_comfixed(); #endif - check_forces(); } @@ -392,27 +391,9 @@ void init_forces() #endif } -MDINLINE void check_particle_force(Particle *part) -{ - - int i; - for (i=0; i< 3; i++) { - if isnan(part->f.f[i]) { - char *errtext = runtime_error(128); - ERROR_SPRINTF(errtext,"{999 force on particle was NAN.} "); - } - } - -#ifdef ROTATION - for (i=0; i< 3; i++) { - if isnan(part->f.torque[i]) { - char *errtext = runtime_error(128); - ERROR_SPRINTF(errtext,"{999 force on particle was NAN.} "); - } - } -#endif -} +// This function is no longer called from force_calc(). +// The check was moved to rescale_fores() to avoid an additional iteration over all particles void check_forces() { Cell *cell; diff --git a/src/forces.h b/src/forces.h index e04f6bbe302..b788905a084 100644 --- a/src/forces.h +++ b/src/forces.h @@ -644,6 +644,27 @@ MDINLINE void add_force(ParticleForce *F_to, ParticleForce *F_add) #endif } +MDINLINE void check_particle_force(Particle *part) +{ + + int i; + for (i=0; i< 3; i++) { + if isnan(part->f.f[i]) { + char *errtext = runtime_error(128); + ERROR_SPRINTF(errtext,"{999 force on particle was NAN.} "); + } + } + +#ifdef ROTATION + for (i=0; i< 3; i++) { + if isnan(part->f.torque[i]) { + char *errtext = runtime_error(128); + ERROR_SPRINTF(errtext,"{999 force on particle was NAN.} "); + } + } +#endif +} + /*@}*/ #endif diff --git a/src/integrate.c b/src/integrate.c index db82131212d..a31cd47a3ea 100644 --- a/src/integrate.c +++ b/src/integrate.c @@ -700,6 +700,7 @@ void rescale_forces() p = cell->part; np = cell->n; for(i = 0; i < np; i++) { + check_particle_force(&p[i]); p[i].f.f[0] *= scale/PMASS(p[i]); p[i].f.f[1] *= scale/PMASS(p[i]); p[i].f.f[2] *= scale/PMASS(p[i]); @@ -730,6 +731,7 @@ void rescale_forces_propagate_vel() p = cell->part; np = cell->n; for(i = 0; i < np; i++) { + check_particle_force(&p[i]); /* Rescale forces: f_rescaled = 0.5*dt*dt * f_calculated * (1/mass) */ p[i].f.f[0] *= scale/PMASS(p[i]); p[i].f.f[1] *= scale/PMASS(p[i]); From 0ba31b364764bd336fccab72c05d5d1cac3cde83 Mon Sep 17 00:00:00 2001 From: Rudolf Weeber Date: Tue, 10 Jan 2012 14:38:11 +0100 Subject: [PATCH 7/8] Bugfix to thermostat: Fixed particles got thermostat forces --- src/thermostat.h | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/thermostat.h b/src/thermostat.h index 67d92538ca4..288cc8e3d48 100644 --- a/src/thermostat.h +++ b/src/thermostat.h @@ -162,8 +162,7 @@ MDINLINE void friction_thermo_langevin(Particle *p) for ( j = 0 ; j < 3 ; j++) { #ifdef EXTERNAL_FORCES -// if (!(p->l.ext_flag & COORD_FIXED(j))) - if (1==1) + if (!(p->l.ext_flag & COORD_FIXED(j))) #endif { #ifdef LANGEVIN_PER_PARTICLE From e9c5466aa8c01b55067f0757d1b2717e2a66e01c Mon Sep 17 00:00:00 2001 From: Olaf Lenz Date: Tue, 10 Jan 2012 16:41:07 +0100 Subject: [PATCH 8/8] Fixed ROTATION_PER_PARTICLE. --- NEWS | 6 ++++ src/communication.c | 4 +-- src/config.h | 1 - src/energy.h | 2 +- src/magnetic_non_p3m_methods.c | 65 +++++++++++++--------------------- src/particle_data.c | 14 ++++---- src/particle_data.h | 2 +- src/pressure.c | 6 ++-- src/rotation.c | 11 +++--- 9 files changed, 51 insertions(+), 60 deletions(-) diff --git a/NEWS b/NEWS index 0227c00bdab..3cf23b87b53 100644 --- a/NEWS +++ b/NEWS @@ -8,6 +8,12 @@ ESPResSo 3.1 User-visible changes -------------------- +* Added new feature LANGEVIN_PER_PARTICLE that allows to set tha + Langevon parameters temperature and gamma per particle. + +* Added new feature ROTATION_PER_PARTICLE that allow to choose whether + a particle has rotational degrees of freedom or not. + * Added new constraint and LB boundary condition "rhomboid". * Renamed Coulomb method maggs to MEMD (inter coulomb maggs => inter diff --git a/src/communication.c b/src/communication.c index aeaeb9fd49f..6e2c82ad60e 100644 --- a/src/communication.c +++ b/src/communication.c @@ -955,7 +955,7 @@ void mpi_send_vs_relative_slave(int pnode, int part) void mpi_send_rotation(int pnode, int part, int rot) { -#ifdef SWITCHABLE_ROTATION +#ifdef ROTATION_PER_PARTICLE mpi_issue(REQ_SET_ROTATION, pnode, part); if (pnode == this_node) { @@ -972,7 +972,7 @@ void mpi_send_rotation(int pnode, int part, int rot) void mpi_send_rotation_slave(int pnode, int part) { -#ifdef SWITCHABLE_ROTATION +#ifdef ROTATION_PER_PARTICLE if (pnode == this_node) { Particle *p = local_particles[part]; MPI_Status status; diff --git a/src/config.h b/src/config.h index 1c5cfb4fada..25786bc9420 100644 --- a/src/config.h +++ b/src/config.h @@ -256,7 +256,6 @@ #define M_SQRT1_2 0.7071067811865475244008443621048490L /* 1/sqrt(2) */ #endif ->>>>>>> master /********************************************/ /* \name exported functions of config.c */ /********************************************/ diff --git a/src/energy.h b/src/energy.h index f1751915466..79784b138f7 100644 --- a/src/energy.h +++ b/src/energy.h @@ -431,7 +431,7 @@ MDINLINE void add_kinetic_energy(Particle *p1) energy.data.e[0] += (SQR(p1->m.v[0]) + SQR(p1->m.v[1]) + SQR(p1->m.v[2]))*PMASS(*p1); #ifdef ROTATION -#ifdef SWITCHABLE_ROTATION +#ifdef ROTATION_PER_PARTICLE if (p1->p.rotation) #endif { diff --git a/src/magnetic_non_p3m_methods.c b/src/magnetic_non_p3m_methods.c index ed357fdd0c7..1c044188536 100644 --- a/src/magnetic_non_p3m_methods.c +++ b/src/magnetic_non_p3m_methods.c @@ -34,7 +34,6 @@ #include "domain_decomposition.h" #include "magnetic_non_p3m_methods.h" - #ifdef DIPOLES // Calculates dipolar energy and/or force between two particles @@ -134,7 +133,18 @@ int tclcommand_inter_magnetic_parse_dawaanr(Tcl_Interp * interp, int argc, char if (coulomb.Dmethod != DIPOLAR_ALL_WITH_ALL_AND_NO_REPLICA ) { coulomb.Dmethod = DIPOLAR_ALL_WITH_ALL_AND_NO_REPLICA; } - + +#ifdef PARTIAL_PERIODIC + if(PERIODIC(0) == 0 || + PERIODIC(1) == 0 || + PERIODIC(2) == 0) + { + Tcl_AppendResult(interp, "Need periodicity (1,1,1) with DAWAANR", + (char *) NULL); + return TCL_ERROR; + } +#endif + if (n_nodes > 1) { Tcl_AppendResult(interp, "sorry: DAWAANR only works with 1 cpu", (char *) NULL); return TCL_ERROR; @@ -152,50 +162,11 @@ double dawaanr_calculations(int force_flag, int energy_flag) double u; int i,j,c,cc; - return 0.0; -} - -/************************************************************/ - - -double dawaanr_calculations(int force_flag, int energy_flag) { - Cell *cell; - Particle *part; - int i,c,np; - double *x=NULL, *y=NULL, *z=NULL; - double *mx=NULL, *my=NULL, *mz=NULL; - double *fx=NULL, *fy=NULL, *fz=NULL; -#ifdef ROTATION - double *tx=NULL, *ty=NULL, *tz=NULL; -#endif - int dip_particles,dip_particles2,j; - double u,r,pe1,pe2,pe3,r3,r5,r2,r7,rx,ry,rz,a,b,cc,d; -#ifdef ROTATION - double bx,by,bz,ax,ay,az; -#endif - double ffx,ffy,ffz; - if(n_nodes!=1) {fprintf(stderr,"error: DAWAANR is just for one cpu .... \n"); exit(1);} if(!(force_flag) && !(energy_flag) ) {fprintf(stderr," I don't know why you call dawaanr_caclulations with all flags zero \n"); return 0;} // Variable to sum up the energy u=0; - for(i=0;ip.rinertia[1] = 1.0; part->p.rinertia[2] = 1.0; #endif -#ifdef SWITCHABLE_ROTATION +#ifdef ROTATION_PER_PARTICLE part->p.rotation =0; #endif @@ -560,7 +560,7 @@ void tclcommand_part_print_virtual(Particle *part, char *buffer, Tcl_Interp *int } #endif -#ifdef SWITCHABLE_ROTATION +#ifdef ROTATION_PER_PARTICLE void tclcommand_part_print_rotation(Particle *part, char *buffer, Tcl_Interp *interp) { sprintf(buffer,"%i", part->p.rotation); @@ -846,7 +846,7 @@ int tclprint_to_result_Particle(Tcl_Interp *interp, int part_num) tclcommand_part_print_torque(&part, buffer, interp); #endif -#ifdef SWITCHABLE_ROTATION +#ifdef ROTATION_PER_PARTICLE Tcl_AppendResult(interp, " rotation ", (char *)NULL); tclcommand_part_print_rotation(&part, buffer, interp); #endif @@ -1008,7 +1008,7 @@ int tclcommand_part_parse_print(Tcl_Interp *interp, int argc, char **argv, else if (ARG0_IS_S("tbf")) tclcommand_part_print_torque_body_frame(&part, buffer, interp); #endif -#ifdef SWITCHABLE_ROTATION +#ifdef ROTATION_PER_PARTICLE else if (ARG0_IS_S("rotation")) tclcommand_part_print_rotation(&part, buffer, interp); #endif @@ -1183,7 +1183,7 @@ int tclcommand_part_parse_rotational_inertia(Tcl_Interp *interp, int argc, char } #endif -#ifdef SWITCHABLE_ROTATION +#ifdef ROTATION_PER_PARTICLE int tclcommand_part_parse_rotation(Tcl_Interp *interp, int argc, char **argv, int part_num, int * change) { @@ -2097,7 +2097,7 @@ int tclcommand_part_parse_cmd(Tcl_Interp *interp, int argc, char **argv, err = tclcommand_part_parse_rotational_inertia(interp, argc-1, argv+1, part_num, &change); #endif -#ifdef SWITCHABLE_ROTATION +#ifdef ROTATION_PER_PARTICLE else if (ARG0_IS_S("rotation")) err = tclcommand_part_parse_rotation(interp, argc-1, argv+1, part_num, &change); #endif @@ -2378,7 +2378,7 @@ int set_particle_rotational_inertia(int part, double rinertia[3]) #endif -#ifdef SWITCHABLE_ROTATION +#ifdef ROTATION_PER_PARTICLE int set_particle_rotation(int part, int rot) { int pnode; diff --git a/src/particle_data.h b/src/particle_data.h index e8144a3c72f..ec1b23d9999 100644 --- a/src/particle_data.h +++ b/src/particle_data.h @@ -79,7 +79,7 @@ typedef struct { double rinertia[3]; #endif -#ifdef SWITCHABLE_ROTATION +#ifdef ROTATION_PER_PARTICLE // Determines, wether a particle's rotational degrees of freedom are // integrated int rotation; diff --git a/src/pressure.c b/src/pressure.c index ce323c22933..35019164c34 100644 --- a/src/pressure.c +++ b/src/pressure.c @@ -139,6 +139,7 @@ void pressure_calc(double *result, double *result_t, double *result_nb, double * virials.data.e[0] /= (6.0*volume*time_step*time_step); #else virials.data.e[0] /= (3.0*volume*time_step*time_step); +#endif calc_long_range_virials(); @@ -147,8 +148,9 @@ void pressure_calc(double *result, double *result_t, double *result_nb, double * /* stress tensor part */ - /* ROTATION option does not effect stress tensor calculations since rotational - energy is not included in the ideal term (unlike for the pressure) */ + /* The ROTATION option does not effect stress tensor calculations + since rotational energy is not included in the ideal term (unlike + for the pressure) */ for(i=0; i<9; i++) p_tensor.data.e[i] /= (volume*time_step*time_step); diff --git a/src/rotation.c b/src/rotation.c index 4d78592be6f..ea301b3df5b 100644 --- a/src/rotation.c +++ b/src/rotation.c @@ -1,6 +1,7 @@ /* - Copyright (C) 2010 The ESPResSo project - Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010 Max-Planck-Institute for Polymer Research, Theory Group, PO Box 3148, 55021 Mainz, Germany + Copyright (C) 2010,2011 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. @@ -251,7 +252,7 @@ void propagate_omega_quat() p = cell->part; np = cell->n; for(i = 0; i < np; i++) { -#ifdef SWITCHABLE_ROTATION +#ifdef ROTATION_PER_PARTICLE if (!p[i].p.rotation) continue; #endif @@ -304,7 +305,7 @@ void convert_torques_propagate_omega() p = cell->part; np = cell->n; for(i = 0; i < np; i++) { -#ifdef SWITCHABLE_ROTATION +#ifdef ROTATION_PER_PARTICLE if (!p[i].p.rotation) continue; #endif @@ -385,7 +386,7 @@ void convert_initial_torques() p = cell->part; np = cell->n; for(i = 0; i < np; i++) { -#ifdef SWITCHABLE_ROTATION +#ifdef ROTATION_PER_PARTICLE if (!p[i].p.rotation) continue; #endif