From debb39968cef240dec8e29bb9232e6c165b0f677 Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Mon, 4 Jul 2011 16:48:18 +0200 Subject: [PATCH 1/4] possibility to keep longrange part of forces and random forces from the thermostat. --- src/forces.c | 18 +++++++++++++++++- src/p3m-magnetostatics.c | 18 ++++++++++++++++++ src/p3m.c | 3 +++ src/particle_data.h | 16 +++++++++++++++- 4 files changed, 53 insertions(+), 2 deletions(-) diff --git a/src/forces.c b/src/forces.c index 1f32982a6cb..c0416782a25 100644 --- a/src/forces.c +++ b/src/forces.c @@ -264,12 +264,23 @@ MDINLINE void init_local_particle_force(Particle *part) part->p.adress_weight=new_weight; } #endif - if ( thermo_switch & THERMO_LANGEVIN ) + if ( thermo_switch & THERMO_LANGEVIN ){ friction_thermo_langevin(part); + #ifdef KEEP_RANDOM + part->l.f_rand[0] = part->f.f[0]; + part->l.f_rand[1] = part->f.f[1]; + part->l.f_rand[2] = part->f.f[2]; + #endif + } else { part->f.f[0] = 0; part->f.f[1] = 0; part->f.f[2] = 0; + #ifdef KEEP_KPART + part->l.f_k[0] = 0; + part->l.f_k[1] = 0; + part->l.f_k[2] = 0; + #endif } #ifdef EXTERNAL_FORCES @@ -287,6 +298,11 @@ MDINLINE void init_local_particle_force(Particle *part) part->f.torque[0] = 0; part->f.torque[1] = 0; part->f.torque[2] = 0; + #ifdef KEEP_KPART + part->l.torque_k[0] = 0; + part->l.torque_k[1] = 0; + part->l.torque_k[2] = 0; + #endif /* and rescale quaternion, so it is exactly of unit length */ scale = sqrt( SQR(part->r.quat[0]) + SQR(part->r.quat[1]) + diff --git a/src/p3m-magnetostatics.c b/src/p3m-magnetostatics.c index d65bc3b51e6..393684a0ebd 100644 --- a/src/p3m-magnetostatics.c +++ b/src/p3m-magnetostatics.c @@ -1125,14 +1125,26 @@ Since the torque is the dipole moment cross-product with E, we have: case 0: //E_x p[i].f.torque[1] -= p[i].r.dip[2]*prefac*dp3m.ca_frac[cf_cnt]*dp3m.rs_mesh[q_ind]; p[i].f.torque[2] += p[i].r.dip[1]*prefac*dp3m.ca_frac[cf_cnt]*dp3m.rs_mesh[q_ind]; + #ifdef KEEP_KPART + p[i].l.torque_k[1] -= p[i].r.dip[2]*prefac*dp3m.ca_frac[cf_cnt]*dp3m.rs_mesh[q_ind]; + p[i].l.torque_k[2] += p[i].r.dip[1]*prefac*dp3m.ca_frac[cf_cnt]*dp3m.rs_mesh[q_ind]; + #endif break; case 1: //E_y p[i].f.torque[0] += p[i].r.dip[2]*prefac*dp3m.ca_frac[cf_cnt]*dp3m.rs_mesh[q_ind]; p[i].f.torque[2] -= p[i].r.dip[0]*prefac*dp3m.ca_frac[cf_cnt]*dp3m.rs_mesh[q_ind]; + #ifdef KEEP_KPART + p[i].l.torque_k[0] += p[i].r.dip[2]*prefac*dp3m.ca_frac[cf_cnt]*dp3m.rs_mesh[q_ind]; + p[i].l.torque_k[2] -= p[i].r.dip[0]*prefac*dp3m.ca_frac[cf_cnt]*dp3m.rs_mesh[q_ind]; + #endif break; case 2: //E_z p[i].f.torque[0] -= p[i].r.dip[1]*prefac*dp3m.ca_frac[cf_cnt]*dp3m.rs_mesh[q_ind]; p[i].f.torque[1] += p[i].r.dip[0]*prefac*dp3m.ca_frac[cf_cnt]*dp3m.rs_mesh[q_ind]; + #ifdef KEEP_KPART + p[i].l.torque_k[0] -= p[i].r.dip[1]*prefac*dp3m.ca_frac[cf_cnt]*dp3m.rs_mesh[q_ind]; + p[i].l.torque_k[1] += p[i].r.dip[0]*prefac*dp3m.ca_frac[cf_cnt]*dp3m.rs_mesh[q_ind]; + #endif } q_ind++; cf_cnt++; @@ -1183,6 +1195,12 @@ static void dp3m_assign_forces_dip(double prefac, int d_rs) ( dp3m.rs_mesh_dip[0][q_ind]*p[i].r.dip[0] +dp3m.rs_mesh_dip[1][q_ind]*p[i].r.dip[1] +dp3m.rs_mesh_dip[2][q_ind]*p[i].r.dip[2]); + #ifdef KEEP_KPART + p[i].l.f_k[d_rs] += prefac*dp3m.ca_frac[cf_cnt]* + ( dp3m.rs_mesh_dip[0][q_ind]*p[i].r.dip[0] + +dp3m.rs_mesh_dip[1][q_ind]*p[i].r.dip[1] + +dp3m.rs_mesh_dip[2][q_ind]*p[i].r.dip[2]); + #endif q_ind++; cf_cnt++; } diff --git a/src/p3m.c b/src/p3m.c index cdbc2ab2df5..b969a359b6f 100644 --- a/src/p3m.c +++ b/src/p3m.c @@ -910,6 +910,9 @@ static void P3M_assign_forces(double force_prefac, int d_rs) for(i1=0; i1 Date: Mon, 4 Jul 2011 17:10:10 +0200 Subject: [PATCH 2/4] added tclcommands and extended part print to output kspace part of forces and torques. --- src/particle_data.c | 42 ++++++++++++++++++++++++++++++++++++++++++ src/rotation.c | 3 +-- src/rotation.h | 4 ++++ 3 files changed, 47 insertions(+), 2 deletions(-) diff --git a/src/particle_data.c b/src/particle_data.c index caec1cdb190..113bbcf5355 100644 --- a/src/particle_data.c +++ b/src/particle_data.c @@ -483,6 +483,27 @@ void tclcommand_part_print_torque(Particle *part, char *buffer, Tcl_Interp *inte Tcl_AppendResult(interp, buffer, (char *)NULL); } +#ifdef KEEP_KPART +void tclcommand_part_print_torque_k(Particle *part, char *buffer, Tcl_Interp *interp) +{ + double torque[3]; +//in Espresso torques are in body-fixed frames. We should convert they to the space-fixed coordinates. + double A[9]; + define_rotation_matrix(part, A); + + torque[0] = A[0 + 3*0]*part->l.torque_k[0] + A[1 + 3*0]*part->l.torque_k[1] + A[2 + 3*0]*part->l.torque_k[2]; + torque[1] = A[0 + 3*1]*part->l.torque_k[0] + A[1 + 3*1]*part->l.torque_k[1] + A[2 + 3*1]*part->l.torque_k[2]; + torque[2] = A[0 + 3*2]*part->l.torque_k[0] + A[1 + 3*2]*part->l.torque_k[1] + A[2 + 3*2]*part->l.torque_k[2]; + + Tcl_PrintDouble(interp, torque[0], buffer); + Tcl_AppendResult(interp, buffer, " ", (char *)NULL); + Tcl_PrintDouble(interp, torque[1], buffer); + Tcl_AppendResult(interp, buffer, " ", (char *)NULL); + Tcl_PrintDouble(interp, torque[2], buffer); + Tcl_AppendResult(interp, buffer, (char *)NULL); +} +#endif + /* tclcommand_part_print_torque_body_frame: function to have the possiblility of printing also the torques in the body_frame to make compatible the manual recovery of saved configurations with the use of the instruction @@ -586,6 +607,19 @@ void tclcommand_part_print_f(Particle *part, char *buffer, Tcl_Interp *interp) Tcl_AppendResult(interp, buffer, (char *)NULL); } +#ifdef KEEP_KSPACE +void tclcommand_part_print_f_k(Particle *part, char *buffer, Tcl_Interp *interp) +{ + /* unscale forces ! */ + Tcl_PrintDouble(interp, part->l.f_k[0]/(0.5*time_step*time_step), buffer); + Tcl_AppendResult(interp, buffer, " ", (char *)NULL); + Tcl_PrintDouble(interp, part->l.f_k[1]/(0.5*time_step*time_step), buffer); + Tcl_AppendResult(interp, buffer, " ", (char *)NULL); + Tcl_PrintDouble(interp, part->l.f_k[2]/(0.5*time_step*time_step), buffer); + Tcl_AppendResult(interp, buffer, (char *)NULL); +} +#endif + void tclcommand_part_print_position(Particle *part, char *buffer, Tcl_Interp *interp) { double ppos[3]; @@ -941,6 +975,10 @@ int tclcommand_part_parse_print(Tcl_Interp *interp, int argc, char **argv, tclcommand_part_print_position(&part, buffer, interp); else if (ARG0_IS_S("force")) tclcommand_part_print_f(&part, buffer, interp); +#ifdef KEEP_KSPACE + else if (ARG0_IS_S("force_k")) + tclcommand_part_print_f_k(&part, buffer, interp); +#endif else if (ARG0_IS_S("folded_position")) tclcommand_part_print_folded_position(&part, buffer, interp); else if (ARG0_IS_S("type")) { @@ -981,6 +1019,10 @@ int tclcommand_part_parse_print(Tcl_Interp *interp, int argc, char **argv, tclcommand_part_print_omega(&part, buffer, interp); else if (ARG0_IS_S("torque")) tclcommand_part_print_torque(&part, buffer, interp); +#ifdef KEEP_KSPACE + else if (ARG0_IS_S("torque")) + tclcommand_part_print_torque_k(&part, buffer, interp); +#endif else if (ARG0_IS_S("tbf")) tclcommand_part_print_torque_body_frame(&part, buffer, interp); #endif diff --git a/src/rotation.c b/src/rotation.c index c0e73f082b9..ca51032fb1f 100644 --- a/src/rotation.c +++ b/src/rotation.c @@ -63,8 +63,7 @@ static double I[3] = { 1, 1, 1}; /************************************************************/ /*@{*/ -/** define rotation matrix A for a given particle */ -static void define_rotation_matrix(Particle *p, double A[9]); + /** define first and second time derivatives of a quaternion, as well as the angular acceleration. */ diff --git a/src/rotation.h b/src/rotation.h index 87e1104eed1..a8067c10ca2 100644 --- a/src/rotation.h +++ b/src/rotation.h @@ -66,4 +66,8 @@ void multiply_quaternions(double a[4], double b[4], double result[4]); int convert_quatu_to_quat(double d[3], double quat[4]); void convert_omega_body_to_space(Particle *p, double *omega); + +/** define rotation matrix A for a given particle */ +void define_rotation_matrix(Particle *p, double A[9]); + #endif From 80ea97b5696efa2cce48d27c9f90bf8531b1d216 Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Mon, 4 Jul 2011 17:22:39 +0200 Subject: [PATCH 3/4] small fix. --- src/particle_data.c | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/particle_data.c b/src/particle_data.c index 113bbcf5355..5e14e1c5138 100644 --- a/src/particle_data.c +++ b/src/particle_data.c @@ -607,7 +607,7 @@ void tclcommand_part_print_f(Particle *part, char *buffer, Tcl_Interp *interp) Tcl_AppendResult(interp, buffer, (char *)NULL); } -#ifdef KEEP_KSPACE +#ifdef KEEP_KPART void tclcommand_part_print_f_k(Particle *part, char *buffer, Tcl_Interp *interp) { /* unscale forces ! */ @@ -975,7 +975,7 @@ int tclcommand_part_parse_print(Tcl_Interp *interp, int argc, char **argv, tclcommand_part_print_position(&part, buffer, interp); else if (ARG0_IS_S("force")) tclcommand_part_print_f(&part, buffer, interp); -#ifdef KEEP_KSPACE +#ifdef KEEP_KPART else if (ARG0_IS_S("force_k")) tclcommand_part_print_f_k(&part, buffer, interp); #endif @@ -1019,7 +1019,7 @@ int tclcommand_part_parse_print(Tcl_Interp *interp, int argc, char **argv, tclcommand_part_print_omega(&part, buffer, interp); else if (ARG0_IS_S("torque")) tclcommand_part_print_torque(&part, buffer, interp); -#ifdef KEEP_KSPACE +#ifdef KEEP_KPART else if (ARG0_IS_S("torque")) tclcommand_part_print_torque_k(&part, buffer, interp); #endif From a093c622e4c8238b59ea0927fd410dfb430b8cee Mon Sep 17 00:00:00 2001 From: Olaf Lenz Date: Thu, 28 Jul 2011 11:26:52 +0200 Subject: [PATCH 4/4] Removed wrong rescaling. --- src/particle_data.c | 7 +++---- testsuite/p3m_wall.tcl | 4 +++- 2 files changed, 6 insertions(+), 5 deletions(-) diff --git a/src/particle_data.c b/src/particle_data.c index 5e14e1c5138..44d547a0265 100644 --- a/src/particle_data.c +++ b/src/particle_data.c @@ -610,12 +610,11 @@ void tclcommand_part_print_f(Particle *part, char *buffer, Tcl_Interp *interp) #ifdef KEEP_KPART void tclcommand_part_print_f_k(Particle *part, char *buffer, Tcl_Interp *interp) { - /* unscale forces ! */ - Tcl_PrintDouble(interp, part->l.f_k[0]/(0.5*time_step*time_step), buffer); + Tcl_PrintDouble(interp, part->l.f_k[0], buffer); Tcl_AppendResult(interp, buffer, " ", (char *)NULL); - Tcl_PrintDouble(interp, part->l.f_k[1]/(0.5*time_step*time_step), buffer); + Tcl_PrintDouble(interp, part->l.f_k[1], buffer); Tcl_AppendResult(interp, buffer, " ", (char *)NULL); - Tcl_PrintDouble(interp, part->l.f_k[2]/(0.5*time_step*time_step), buffer); + Tcl_PrintDouble(interp, part->l.f_k[2], buffer); Tcl_AppendResult(interp, buffer, (char *)NULL); } #endif diff --git a/testsuite/p3m_wall.tcl b/testsuite/p3m_wall.tcl index 25d8c26ac9e..254655cdd66 100644 --- a/testsuite/p3m_wall.tcl +++ b/testsuite/p3m_wall.tcl @@ -71,7 +71,9 @@ integrate 0 set rmsf 0 for { set i 0 } { $i <= [setmd max_part] } { incr i } { - set resF [part $i pr f] + set resF [part $i print f] + puts [part $i print id force] + puts [part $i print id force_k] set tgtF $F($i) set dx [expr ([lindex $resF 0] - [lindex $tgtF 0])] set dy [expr ([lindex $resF 1] - [lindex $tgtF 1])]