Permalink
Browse files

Merge pull request #45 from rempferg/galilei_transform

Center of mass motion tools compatible with MASS
  • Loading branch information...
2 parents 360e18f + f17675b commit c4da8bec6a0abff3e41a3562c98d654a38f074f4 @olenz olenz committed Mar 16, 2012
Showing with 106 additions and 49 deletions.
  1. +29 −4 doc/ug/aux.tex
  2. +4 −5 doc/ug/setup.tex
  3. +73 −40 scripts/auxiliary.tcl
View
33 doc/ug/aux.tex
@@ -21,12 +21,37 @@ \chapter{Auxilliary commands}
\label{chap:aux}
\todo{Missing commands:
- Probably all from \keyword{scripts/auxiliary.tcl}?
- \keyword{galileiTransformParticles},
- \keyword{system_com_vel}}
+ Probably all from \keyword{scripts/auxiliary.tcl}?}
-\section{Finding particles and bonds}
+\section{Center of mass motion}
+\label{sec:centerofmassmotion}
+
+\subsection{\texttt{system_com}}
+\begin{essyntax}
+ system_com
+\end{essyntax}
+Returns the center of mass of the whole system.
+
+\subsection{\texttt{system_com_vel}}
+\begin{essyntax}
+ system_com_vel
+\end{essyntax}
+Returns the velocity of the center of mass of the whole system.
+\subsection{\texttt{galileiTransformParticles}}
+\begin{essyntax}
+ galileiTransformParticles
+\end{essyntax}
+Substracts the velocity of the center of mass of the whole system from every
+particle's velocity, thereby performing a galilei transform into the reference
+frame of the center of mass of the system.
+
+This is useful for example in combination with the DPD thermostat, since there,
+ a drift in the velocity of the whole system leads to an offset in the reported
+ temperature.
+
+
+\section{Finding particles and bonds}
\subsection{\texttt{countBonds}}
\begin{essyntax}
View
9 doc/ug/setup.tex
@@ -271,11 +271,10 @@ \subsubsection{Thermostat DPD}
temperature jumps during a simulation are. The dpd thermostat does not
act on the system center of mass motion. Therefore, before using dpd,
you have to stop the center of mass motion of your system, which you
-can achieve by using the command \keyword{galileiTransformParticles}.
-\todo{Docs for galileiTransformParticles is missing} This may be
-repeated once in a while for long runs due to round off errors (check
-this with the command \keyword{system_com_vel}). \todo{Docs from
- system_com_vel is missing}
+can achieve by using the command \keyword{galileiTransformParticles}
+\ref{sec:centerofmassmotion}. This may be repeated once in a while for long
+runs due to round off errors (check this with the command
+\keyword{system_com_vel}) \ref{sec:centerofmassmotion}.
Two restrictions apply for the dpd implementation of \es:
View
113 scripts/auxiliary.tcl
@@ -467,51 +467,83 @@ proc stopParticles { } {
# center of mass
# --------------
#
-# Calculate the center of mass motion of the system stored in Espresso
+# Calculate the center of mass of the system stored in Espresso
#
#############################################################
proc system_com { } {
- set npart [setmd n_part]
+ set npart [setmd n_part]
- # calculate center of mass motion
- set com { 0 0 0 }
- set part_cnt 0
- set i 0
- while { $part_cnt < $npart } {
- if { [part $i] != "na" } {
- set pos [part $i print p]
- set com [vecadd $com $pos]
- incr part_cnt
- }
- incr i
+ #calculate center of mass
+ set com {0 0 0}
+ set part_cnt 0
+
+ if {[has_feature MASS]} {
+ set net_mass 0
+
+ for {set i 0} {$part_cnt < $npart} {incr i} {
+ if {[part $i] != "na"} {
+ set pos [part $i print p]
+ set mass [part $i print mass]
+ set com [vecadd $com [vecscale $mass $pos]]
+ set net_mass [expr $net_mass + $mass]
+ incr part_cnt
+ }
+ }
+
+ return [vecscale [expr 1.0/$net_mass] $com]
+ } else {
+ for {set i 0} {$part_cnt < $npart} {incr i} {
+ if {[part $i] != "na"} {
+ set pos [part $i print p]
+ set com [vecadd $com $pos]
+ incr part_cnt
+ }
}
+
return [vecscale [expr 1.0/$npart] $com]
+ }
}
#
# center of mass motion
# ---------------------
#
-# Calculate the center of mass motion of the system stored in Espresso
+# Calculate the center of mass velocity of the system stored in Espresso
#
#############################################################
-proc system_com_vel { } {
- set npart [setmd n_part]
+proc system_com_vel {} {
+ set npart [setmd n_part]
- # calculate center of mass motion
- set com_vel { 0 0 0 }
- set part_cnt 0
- set i 0
- while { $part_cnt < $npart } {
- if { [part $i] != "na" } {
- set part_vel [part $i print v]
- set com_vel [vecadd $com_vel $part_vel]
- incr part_cnt
- }
- incr i
+ # calculate center of mass motion
+ set com_vel {0 0 0}
+ set part_cnt 0
+
+ if {[has_feature MASS]} {
+ set net_mass 0
+
+ for {set i 0} {$part_cnt < $npart } {incr i} {
+ if { [part $i] != "na" } {
+ set part_vel [part $i print v]
+ set part_mass [part $i print mass]
+ set com_vel [vecadd $com_vel [vecscale $part_mass $part_vel]]
+ set net_mass [expr $net_mass + $part_mass]
+ incr part_cnt
+ }
+ }
+
+ return [vecscale [expr 1.0/$net_mass] $com_vel]
+ } else {
+ for {set i 0} {$part_cnt < $npart } {incr i} {
+ if { [part $i] != "na" } {
+ set part_vel [part $i print v]
+ set com_vel [vecadd $com_vel $part_vel]
+ incr part_cnt
+ }
}
+
return [vecscale [expr 1.0/$npart] $com_vel]
+ }
}
#
@@ -521,24 +553,25 @@ proc system_com_vel { } {
# Perform a Galilei Transformation on all Particles stored
# in Espresso such that the center of mass motion of the
# system is zero. Returns old center of mass motion.
-# Assumes particles of equal mass!
#
#############################################################
-proc galileiTransformParticles { } {
- set com_vel [system_com_vel]
+proc galileiTransformParticles {} {
+ set com_vel [system_com_vel]
- # subtract center of mass motion from particle velocities
- set npart [setmd n_part]; set part_cnt 0; set i 0
- while { $part_cnt < $npart } {
- if { [part $i] != "na" } {
- set part_vel [vecsub [part $i print v] $com_vel]
- part $i v [lindex $part_vel 0] [lindex $part_vel 1] [lindex $part_vel 2]
- incr part_cnt
- }
- incr i
+ #subtract center of mass motion from particle velocities
+ set npart [setmd n_part];
+ set part_cnt 0;
+
+ for {set i 0} { $part_cnt < $npart } {incr i} {
+ if {[part $i] != "na"} {
+ set part_vel [vecsub [part $i print v] $com_vel]
+ part $i v [lindex $part_vel 0] [lindex $part_vel 1] [lindex $part_vel 2]
+ incr part_cnt
}
- return $com_vel
+ }
+
+ return $com_vel
}
#

0 comments on commit c4da8be

Please sign in to comment.