Subject: 100% parallel MD(2nd Summary)

From: [email-masked] (Noy)
Date: Tue, 1 Nov 1994 09:14:39 +0100 (NFT)

Dear fellows on the cyberspace,
        This is the summary of the responses I have got. Hope you
enjoy it. Happy Halloween......
                                                       take care,
                                                       Teerakiat
----------------------------  original post -------------------------
Dear cyberchemists,
        From the principle of molecular dynamics, velocity is
reversible. Thus, if we assign a minus sign to all velocity of every
particle in the simulation box ( at any equilibrated state ), the
trajectory should be recoverd to the past correctly. 
        This is not a new idea but I wonder if somebody has
practiced it. For example if we simulate the phase space at time
orgin t0 to the time t1,
        
        t0 ---------------->   t1     (1 picosecond)

        Then, we revert velocity of all particles of the configuration
at the time t0. Simulates it back to the past to the time t2

        t0 ---------------->   t2     (1 picosecond)

        Then we have the past trajectory.
        If we connect the two trajectories together, we will have

         t2 -------------------> t0 -----------------> t1 (2 picoseconds)
       
        Thus, we must revert the sign of velocity of the configurations
>from the time t1 to the time t0. We will have the overall longer trajectory
of the simulation. Given the run of two trajectories running in parallel
on two computers, one can save time by the factor of 2. 
        The dynamics and structural properties should be calculated
>from the time origin t2 to the time t1.
        Does somebody have comments or suggestions on this ? Is it
correct ? If this idea is practical, it would help to save time a lot.
Imagine if I have to run MD for 10 ps to evaluate diffusion coefficient
more precisely, I can divide the work into 2 jobs. One job is simulated
on a computer. For another job, I convert the velocity and trace the
phase space back to it past, have it run on another computer parallelly.
Connect two jobs together by going back to the origin of the phase space
in the far past time. Evaluate the time evolution and then I should
get the results much faster. After all, the parallel efficiency is 
100 %  !!!!!!!!
        Thanks for your comments.
                                                best wishes,
                                                Teerakiat
---------------------------- end of original post -----------------------

---------------------------- responses -----------------------------------
From: Richard A Caldwell <[email-masked]>

Does the reversal risk numerical error? Small but numerous steps have 
roundoff error, and large steps have concomitant algorithmic inaccuracy. 
I am just wondering.

dc

----------
From: [email-masked]

From: Dennis J. Gerson, Ph.D.
      Technical Market Support; Chemical Sciences

The parallel efficiency will NEVER be 100% because of setup time and
shutdown time within the program (Amdahl's Law).  The best you can do is
to get less than 1% serial content in your code since you will need to
have serial code to read in the starting information and write out the
results (Amber is 97% parallel).  Also, you must consider how you plan
to partition the data and how the processes will communicate with
each other.  If the "scatter data and processes" and "gather data and
processes" (either via fork/join on an SMP or send/receive on a DMP
system) require more bandwidth than the communication bus can handle,
then you reach the practical limit of parallelism.

Remember, the algorithm may be 100% parallel in content, but the setup
and shutdown are not, so the overall code is NOT 100% parallel.  On top
of this, the hardware/operating_system_software platform introduce
constraints on parallel efficiency via bandwidth and compiler capability.

Regards, Dennis Gerson

---------------------------------------------------------------
IBM POWER Parallel Systems Division |Tele:(214)406-7267 Fax:7246
1503 LBJ Freeway MS/280750          |IBMLINK: USIB5MSK @ IBMMAIL
Dallas TX 75234  USA                |Email: [email-masked]
*** Forwarding note from SMTP2   --IINUS1   10/27/94 09:15 ***
=========================================================================

----------
From: "Gerald Loeffler" <[email-masked]>

I see basically 2 problems with your suggestion:
	1) Certain algorithmic tricks, like
		o cutting off long-range interactions
		  beyond a certain radius
		o simulating with velocity scaling to keep the temperature
		  constant
	   are said to 'destroy reversibility'. As I understand it, this means
	   that you can't go from state S1 to state S2, invert the velocities
	   at state S2, go back the same time to state S1' and expect S1 = S1'.

	   I'm really not sure if this spoils your idea: Phase space sampling
	   is certainly not affected by this argument and should be perfectly
	   compatible with your suggestion even when using cut-offs and
	   velocity-scaling. Whether time correlations have an discontinuity
	   at the point where you start simulating in two directions is
	   not obvious.

	2) The system has to be equilibrated, which it usually isn't when
	   you start a MD. But this just means that you would have to simulate
	   - say - 100ps as usual, and then apply your method.

I think that the idea is very interesting, and would like to hear what other
answers you got!

	Yours,
	Gerald
--
Gerald Loeffler : [email-masked]

Institute for Theoretical Chemistry
Theoretical Biochemistry Group
University of Vienna
Waehringerstrasse 17/Parterre
A-1090 Wien, Austria

----------

From: DiveMaster <[email-masked]>

Suffice it to say, I believe that this is the concept employed 
in parallel processing computers, although your particular 
application of "job-splitting" is equation (whatever you are 
doing) dependent.

				Brent Bettini
				NC School of Science & Math
				Durham, NC  27715
				[email-masked]
----------
From: [email-masked] (David Case)

> 
>         t0 ---------------->   t2     (1 picosecond)
> 
>         Then we have the past trajectory.
>         If we connect the two trajectories together, we will have
> 
>          t2 -------------------> t0 -----------------> t1 (2 picoseconds)
> 
Note that, because of round-off and integration errors (arising from
a finite time step) with most integrators if you actually started at
time t2, you would not reproduce the above trajectory.  But that's
probably not real important in most cases.
       
 After all, the parallel efficiency is 
> 100 %  !!!!!!!!

The idea is sound, but perhaps not worth eight explanation points :-).
You do have to initialize two MD runs rather than one, so the actual
efficiency is less than 100%.  Plus you have to have a way to reverse
the record of one of the trajectories, which can be cumbersome for
long trajectories.  But the main problem is that this "parallel" code
only works on exactly two nodes.

...dave case
 
----------

From: [email-masked]

Dear Teerakiat

	I think that you are forgetting something in your scheme. The fact that
the velocity is a function of the coordinates. In this way, these two steps
that share the same absolute velocity, are only valid at one integration time
step. After that, the coordinates of the system are no longer the same.

	Regards,
	Claudio
-------------------------------------
Phone : (351-1) 442-6146 Ext. 237/341
  Fax : (351-1) 442-8766
email : [email-masked] 
Mailing address :
Claudio Soares
Instituto de Tecnologia Quimica e Biologica
Rua da Quinta Grande 6
Apartado 127
2780 OEIRAS Portugal

----------

From: [email-masked] (James A Given)

  MD trajectories are not reversible if they are longer than a few
time steps. The reason is that the trajectories have nontrivial
"mixing" properties and thus they lose information about the
past. A theorist named Ronald Fox at Georgia Institute of
Technology studied this question and wrote about it in Physical 
Review A (late 1970's ??). Do the experiment. It doesn't
take very long. Run a trajectory. Reverse it. See what results.

----------

From: [email-masked] (Ferenc Molnar)

Dear Teerakiat:
Even if your integrator provides equations which are reversible in time, you
should always worry about numerical precision. Try to check out if the
integration from t1 with all velocities reversed leads to t0 and the
same for t2 to t0. Otherwise your combined trajectory won't be correct.
If this very nice idea turns out to be fruitfull, please let me know.

Cheers,
Ferenc
Ferenc Molnar
---------------------------------------------------------------------------
Institut fuer Physikalische und Theoretische Chemie
- Lehrstuhl Prof. Dick -                   Tel.:  (+49) 941 943-4466 /-4486
Universitaet Regensburg                    Fax.:  (+49) 941 943-4488
Universitaetsstrasse 31
D-93053 Regensburg
Deutschland / Germany
---------------------------------------------------------------------------
EMail (SMTP):  [email-masked]
---------------------------------------------------------------------------
:-) ...THE PHYSICISTS HAVE MADE THEIR UNIVERSE, AND IF YOU DO NOT LIKE
 IT, YOU MUST MAKE YOUR OWN.
                               -- E. R. HARRISON IN "COSMOLOGY" (1980)
---------------------------------------------------------------------------

From: [email-masked]

Dear Klaus

although I am quite new to the field of MD, I would like to share my
thoughts on your posting with you.

It is true that Newton's equation of motion are symmetric with 
respect to time, as they only contain second derivatives. However,
these equations are most often modified in MD calculations.
If one simulates a microcanonical ensemble (constant N,V and total
energy of the system), there is no problem with your approach. 
But, and this is the more practical case, if you simulate a canonical
ensemble (constant N,V and Temperature), you have to make sure that
the temperature remains constant. Often, this is done by weak coupling
of the system to a temperature bath.
This introduces disspative terms into the equations of motion, which makes
the reversal of the time direction impossible.

Besides, if you are interested in longer trajectories to get better
statistics of the properties you derive from your simulation, then you 
might alternatively start a second calculation from a different initial
condition on a second machine.
 
Kind regards,

Marcel Utz.

------------------------------------------------------------------------------
Marcel Utz                                              phone:  +41 1 632 5672 
Institute of Polymers                                     fax:  +41 1 632 1096
ETH-Zurich CNB E 98.2
CH-8092 Zurich, Switzerland                  email: [email-masked]
------------------------------------------------------------------------------

From: [email-masked] (Walter E. Reiher III)

Noy:

You wrote:
>...
>         If we connect the two trajectories together, we will have
> 
>          t2 -------------------> t0 -----------------> t1 (2 picoseconds)
>        
>         Thus, we must revert the sign of velocity of the configurations
> from the time t1 to the time t0. We will have the overall longer trajectory
> of the simulation. Given the run of two trajectories running in parallel
> on two computers, one can save time by the factor of 2. 
>         The dynamics and structural properties should be calculated
> from the time origin t2 to the time t1.
>
>         Does somebody have comments or suggestions on this ? Is it
> correct ? If this idea is practical, it would help to save time a lot.
> Imagine if I have to run MD for 10 ps to evaluate diffusion coefficient
>...
Of course, you have to be certain that 't0' is NOT the beginning of your
MD run, but the beginning of your "collection time", AFTER
equilibration.  Because of integration errors and rescaling of velocity
for temperature equilibration, 'ts' (start of equibration) != 't2':
           ts -------------------> t0
                                  then
           t2 <------------------- t0
                                   t0 -----------------> t1
                            combine to give
           t2 -------------------> t0 -----------------> t1 (2 picoseconds)

The idea makes sense to me, and I've heard people talk about doing this
for years now, but I'm not sure if it's actually been done.  This
technique has, of course, been used to check the integration in MD runs:
after equilibration, it should be possible to run 'forward' then
'backward' in time and get to where you began the 'collection' part (t0
== t2).
           ts -------------------> t0
                                   t0 -----------------> t1
                                   t2 <----------------- t1

For your application, a diffusion coefficient, it sounds reasonable to
me; for other applications of MD, like conformational sampling/
equilibrium structure, it's better to use multiple starting structures
and multiple parallel runs (where 'multiple' can be > 2).

Good luck!

Wally
========================================================================
Walter E. Reiher III, Ph.D.                            [email-masked]
Consultant in Computational Chemistry
P.O. Box 61056                                        voice 408-720-0240
Sunnyvale, CA 94088                                     fax 408-720-0378

---------------
From: K Bryson <[email-masked]>

	Hi Teerakiat,

		Took a while to work out what you were suggesting,
	but finally got it, and it's a very clever idea for parallising
	into two a single run.

		You would of course have to start with a state which
	had been equilibrated before hand, and this process could only
	be done in serial.
		After initial equilibration you would also have to
	perturb something or otherwise the backward time trajectory
	would just go back to your original starting ( unequilibrium )
	state, well for classical dynamics anyway ( ignoring the effect
	of temperature adjustment, this may actually provide the 
	perturbation that will make it produce a good equilibrium run
	backwards.)

	If the system was coupled to a heat bath then it wouldn't go
	back to the starting state (unequilibrated ), but this brings
	in the question whether the temperature coupling process is still
	valid for 'backward' trajectories in time ... certaining if it
	was temp. coupled your splitting process would not end up
	producing the same trajectory as a single 2t length trajectory
	starting from point t2, I think.

		So I think it's a good suggesting if you have an
	equilibrated state and you are running classical dynamics.
	Do not think it will work for temperature coupled systems,
	well, it wouldn't give you an identical path to the 2t length
	simulation, although it might be just as physically realistic
	for calculating time statistics data, not sure.

		What do you think ?

		Interesting idea though,
					Kev.

                              www
                             (o o)
=======================---ooO-(_)-Ooo---=====================================
                               
    K.Bryson                 email: [email-masked]   
    Biophysics Group         Tel  : +44 904 430000 Extn. 2236  
    Physics Department       Fax  : +44 904 432214
    University of York      
    Heslington              "Molecular modelling of DNA and its
    YORK, UK                    interaction with small molecules."
    YO1 5DD             

=============================================================================

From: [email-masked] (Song Ling)

You proposed that by reversing all momenta (velocities when they are 
proportional to momenta) you can connect two pieces of integrated 
trajectories together to form a longer one, unfortunately this is 
incorrect in general.  Here is a simple example, the trajectory of 
a harmonic oscillator in phase space is a closed topological circle,
and for an initial condition (x,p) in the first quadron (spell), 
when you reverse its momentum (velocity), you will end up in the 
4th quadron, integrating the trajectory for a time dt, your first 
trajectory is still in the first quadron, while your another trajectory 
is still in the 4th quadron, the two pieces simply don't connect.

The problem originated from that you didn't think about the trajectories
as phase space objects (flows), you had a coordinate space picture.

It is indeed possible to connect two pieces of trajectories to form a 
longer one in phase space, the way to do it is to integrate backward in
time for -dt, and then integrate forward in time for dt, endding up a
trajectory of "length" 2dt.  Do not reverse the signs of the momenta 
(velocities) of your initial point.  (You stay at the same point in 
phase space.)  But this approach may be error prone depending on how 
many time steps you integrate your trajectory, if your system is very
sensitive to initial conditions, the chances are after you integrate a
trajectory from point A to point B, your integrate backward from B
(with -dt) you may not end up in A.

You seemed to be a little confused with the concept of "time reversal
symmetry", to be safe let's say you have a conservative system (non-
dissipative), the total energy is a constant of the motion, time 
reversal symmetry says that looking at the trajectories in phase space
you don't find the arrow of time.  In the harmonic oscillator example
above, you have a topological circel in the phase space, you can 
parameterize it such that a point on it runs counterclockwise or 
clockwise, corresponding to running time forward or backward.  If you
have a damped harmonic oscillator, your trajectory spirals into the 
origin, the fixed point, now you have an arrow of time, past is mapped 
into the future.

In short, you do not apply time reversal to join two pieces of trajectories
together; you integrate forward and backward in time to do that.  I don't
know if that will save you computer time.

-----------------
From: Gustavo Mercier <[email-masked]>

On Thu, 27 Oct 1994, Noy wrote:

...
>         This is not a new idea but I wonder if somebody has
> practiced it. For example if we simulate the phase space at time
> orgin t0 to the time t1,
>         
>         t0 ---------------->   t1     (1 picosecond)
> 
>         Then, we revert velocity of all particles of the configuration
> at the time t0. Simulates it back to the past to the time t2
> 
>         t0 ---------------->   t2     (1 picosecond)
> 
>         Then we have the past trajectory.
>         If we connect the two trajectories together, we will have
> 
>          t2 -------------------> t0 -----------------> t1 (2 picoseconds)
>        
>         Thus, we must revert the sign of velocity of the configurations
> from the time t1 to the time t0. We will have the overall longer trajectory
> of the simulation. Given the run of two trajectories running in parallel
> on two computers, one can save time by the factor of 2. 
>         The dynamics and structural properties should be calculated
> from the time origin t2 to the time t1.

Hi!

I think there is one problem. By following your prescription you fail
to properly sample the phase space. In fact it looks like you sample
the same phase space backwards and forwards, hence your computation
of dynamic and structural properties will be off. The key is not the
length of the simulation, but the sampling of the phase space. If you
do so in one picosecond, a one picosecond trajectory will be adequate
to compute the property of interest.

I recommend an article by van Gunsteren and Berendsen (Angew. Chem. Int. Ed.
Engl. 29, (1990) 992-1023. He argues that your simulations should not
sample those degrees of freedom that are averaged within the time
frame of the property of interest. Instead the effect of such averaged
motions should be included through some potential of mean field.
In this way you don't have to waste simulation time in sampling the
phase space so that the proper averaging occurs. In other words,
don't do md simulations if your interest is molecular brownian motions.

Good Luck!

Gustavo A. Mercier, Jr.
[email-masked]

-------------------
From: Martyn Winn <[email-masked]>

Teerakiat,
My feeling is that your scheme should work, provided
you also reverse quantities such as the 3rd derivative
of the positions (needed in a predictor-corrector algorithm).
It should then be entirely equivalent to running the
simulation backwards in time.
With regard to Song Ling's objection, he is right that
the two trajectories do not connect in phase space after
the transformation p -> -p (to take his harmonic oscillator
example). However, my understanding is that you would make
the transformation -p -> p after the run, thus reconnecting
the trajectories.

Martyn

--
***************** Dr. Martyn Winn **************************************
 \                                                                    /
  \ Institut fuer Theoretische Physik, Technische Universitaet Wien, /
   \  Wiedner Hauptstrasse 8-10, A-1040 Wien, Austria.              /
    \____                                              ____________/
         \      Tel: +43 1 58801 5678                 /  
          \     Fax: +43 1 5867760                   /
           \    E-mail: [email-masked]     /
            \                                      /
             --------------------------------------

--------------
From: [email-masked] (Song Ling)

Dr. Winn's reconnection scheme is definitely correct, and in principle it
is equivalent to integrating backward as I suggested earlier.  One needs
the 3rd derivatives because predictor-corrector is not self-starting
(unless one uses Runge-Kutta which is slow).  Two questions remain, one
is computational: Will the -dt, dt approach be more efficient than the
"conventional" 2dt forward method?  The other is "fundamental": Will the
-dt, dt approach generate the "same" trajectory as the 2dt one if the
system is sensitive to initial conditions?  I'd like to learn about the 
answers.

-----------------
From: Martyn Winn <[email-masked]>

Two points:
1) My impression is that there is much work devoted to
chaotic trajectories for small molecules (i.e. with a
small no. of degrees of freedom) but that the issue is
ignored for MD simulations of bulk systems. This is, of
course, purely for reasons of convenience.
2) If one is in a chaotic regime, then I'm sure that
the -dt, dt trajectory will differ from the 2dt trajectory,
but taken individually each would be a perfectly valid
trajectory. There would be no reason to prefer one over
the other. In fact, a chaotic trajectory is surely ideal
to satisfy the ergodic requirements of most MD simulations.

I wonder if anyone is answering the original question
about efficiency?!

Martyn Winn

--
***************** Dr. Martyn Winn **************************************
 \                                                                    /
  \ Institut fuer Theoretische Physik, Technische Universitaet Wien, /
   \  Wiedner Hauptstrasse 8-10, A-1040 Wien, Austria.              /
    \____                                                          /
         \      Tel: +43 1 58801 5678                             /  
          \     Fax: +43 1 5867760                               /
           \    E-mail: [email-masked]                 /
            \   URL: http://tph.tuwien.ac.at/~winn/me.html     /
             --------------------------------------------------

----------------
From: Gustavo Mercier <[email-masked]>

Hi!

In light of recent messages, I reviewed the message I sent you.
I realized that I made a mistake that invalidates my argument.
You reverse the simulation at t0 to go to t2, not at t1!
Hence, you would sample a different area of the phase space.

Good Luck!

Gustavo A. Mercier, Jr.
[email-masked]
-------------------

From: "Marais, Charles F. "  <[email-masked]>

> Will the
> -dt, dt approach generate the "same" trajectory as the 2dt one if the
> system is sensitive to initial conditions? 

This is very interesting.  Has much work been done on chaos in MD ? 
 Any pointers would be much appreciated.

C

Dr Charles Marais

Department of Chemistry
University of Cape Town
Private Bag
Rondebosch                                  [email-masked]
South Africa 7700                           [email-masked]
------------------------------------------------------------------