ProtoMol Research Log for
Thierry Matthey
email me at: tmatthey@nd.edu

I added spherical boundary conditions that consist of a single
potential functions. It is implemented as a force and the user adds it
to the integrator definition to impose spherical boundary conditions.
 

Integrator {
  level 0 Leapfrog {
        timestep 1
    force Spherical -center 0 0 0 -radius 9 -k 1 -j 2
    ...

  }
}

The energy contribution are added to 'E_other'. The spherical boundary
conditions include a parallel implementation. It's only supported for
vacuum (normal boundary conditions).

Example:
examples/water_423/equil298K_01.conf.sphericalBC

The equation are from:
http://www.ks.uiuc.edu/Research/namd/1.5/ug/node58.html

Description at:
http://www.ks.uiuc.edu/Research/namd/1.5/ug/node38.html

I did some changes in the multi-grid method, such that you have to
specify the mesh size of the finest grid for each dimension in case of
vacuum (normal boundary conditions). Moreover, I made the re-size of
the grid more conservative since the energy does not retain the
property of being invariant under global translation and/or rotation
of particle positions. I ran the 432-water test case for 1ns, dt=1fs
with a total energy drift of less than 0.3%. The diffusion was
app. 100, due to the fact that multi-grid does not preserve linear
momentum. The grid at finest level was moved 76 times from originally
(-15,-18,-18)-(21,18,18) to (12,234,6)-(48,270,42). The maximum
relative force error is app. 3%.

 Coulomb -algorithm MultiGrid -interpolation Hermite -kernel C2
   -levels 2
   -s 6.5
   -order 4
   -ratio 2
   -h 3 3 3
   -origin 0 0 0

By random I found a glitch in the cell enumerator for PBC when (2*cutoff >
cell dimension). It may be slower but it works for arbitrary cutoff.
Additionally Ewald and PME did not handle (2*Rc > cell dimension) correctly.
The lattice vectors (n != 0) where left out.

For vacuum, MG expects the mesh sizes of the finest grid. It figures out the number of grid points for each level. I made the re-size more conservative, since there is a energy drift when re-sizing or moving the gridpoints.

I did some validation test for Ewald and PME with water_423 (25x25x25 AA^3)

For Ewald it is pretty easy to control the accuracy. As you can see the
relative max error corresponds the chosen accuracy. On the IBM with 1e-16
precision you set the accuracy in the reference Ewald force to 1e-18. It's
partially the same as for a run with accuracy 1e-50. If you are still
not convinced, you run with normal boundary conditions and compare
against the direct method:

force compare time Coulomb -algorithm FullEwald -correction -real -reciprocal
      -accuracy 1e-50
      -j 100000000
force compare time Coulomb -algorithm NonbondedSimpleFull

The relative max error is less than 5e-15. That is what you can expect
due to the machine precision.
 

For PME it is a little more compilcated, but I was able to attain a
relative max error less than 2e-14:

force compare time Coulomb -algorithm PMEwald -correction -real -reciprocal
      -interpolation BSpline
      -gridsize 250 250 250    #mesh size 0.1AA
      -cutoff 10
      -order 12
      -accuracy 1e-18
force compare time Coulomb -algorithm FullEwald -correction -real -reciprocal
      -accuracy 1e-18

Note, that these comparisons are with exclusions. The reciprocal term of
PME is not a negligible number (E-PME=693.2774, E-total=-806.9964).
 

I revisited the timer functions. The framework comes with a nice Timer
class and TimerStatistic class such that the timing is decoupled from
the Parallel class -- the way it should be. The timer comes with a
TimeRep class, which is more or less a container to keep real, user
and system time. Timer class implmenets a timer and does also read the
time from the system. Furthermore I added a wrapper to measure the
time for any desired force and I removed this feature in the
comparisons class CompareForce. However, you can still compare two
forces and measure the time just by combining them:

force compare time ...

You find an example in examples/water_423/equil298K_01.conf.MultiGrid
Note that you should first do the timing and then the comparison, when
combining them ... otherwise you will time also the comparison.

For the overall run time (wall) the real, user and system time is
printed, where for the run, integration, forces, communication and
idle the process time is displayed. The process time is the system
time plus the user time. For me it seems more fair to include also the
system time, because it reflects the time on an empty machine, and
system time is also time protomol has spent.

How does the timers work?

// Starts the timer for the wall clock
TimerStatistic::MyTimer[TimerStatistic::WALL].start();

// Stops the timer for the wall clock
TimerStatistic::MyTimer[TimerStatistic::WALL].stop();
 

... or you can create your own timer:

#include "Timer.h"
Timer timer;
timer.start();
//
// Here I execute my stuff I want to measure
//
timer.stop();
report << plain <<"Time:"<<timer.getTime()<<endr;

If you need the current time you may call this directly by
Timer::getCurrentTime(). It will return a TimeRep object, repesenting
the time (real, user and system).

...Comments...

Performance Test:fire.ii.ui.no IBM/Linux Cluster

  ProtoMol   Namd    Test Case
------------------------------------------------
  10.8163    12.15   1000/NBC/water_423
  15.5792    14.88   1000/PBC/water_423
   7.8260     8.89   10/NBC/bptiAndWater_14281
  13.8554    12.16   10/PBC/bptiAndWater_14281
  56.2211    65.37   10/NBC/apoa1AndWater_92224
  92.5859   100.56   10/PBC/apoa1AndWater_92224

I did some minor changes around the cell manager to improve the
performance, but the real boost you will only get when choosing the right
cell size. For NBC you may try a cell size of 1/1, 1/2+epsilon and
1/3+epsilon of the cutoff. I would recommend 1/2. It may also depend on
your cache and cutoff size. In case of PBC is not really clear what is
best, because the cell size should devied the dimensions of the simulation
box. I would recommend 1/3+epsilon. Both the dimensions  of the
simulation box and the cutoff should be equal n*cutoff-epsilon, where n is
a integer and epsilon as small pos. number. Note that if the cell size is
equal n*cutoff+epsilon you will end with a big overhead due the most outer
cells have only a few pair inside the cutoff. The same yields for the
dimension of the simulation box for PBC. The subdivision into cells starts
in the min corner of the simulation box such that subdivision matches the
boundaries. In the max corner the subdivision will not fit with the
boundaries (only if the cell size is a fraction of the
simulation dimensions). These cells will be added to their closest
neighbors to make the cell enumerator simpler.

Note that pairwise potentials can not relay on i<j any more.

For PME and MG using RefArray should give some boost
inside nested loops. It's already done for the cell manager.

Performance Test:fire.ii.ui.no IBM/Linux Cluster

ProtoMol:

Integrator {
  level 0 Leapfrog {
        timestep 1
    force Improper
    force Dihedral
    force Bond
    force Angle
    force LennardJones Coulomb
        -switchingFunction C2
        -switchingFunction Shift
        -algorithm NonbondedCutoff
        -switchon 0.1
        -cutoff 10.0
  }
}
 

  ProtoMol   Namd    Test Case
------------------------------------------------
  11.2171    12.08   1000/NBC/water_423
  14.6820    14.81   1000/PBC/water_423
  10.5377     9.8    10/NBC/bptiAndWater_14281
  20.6035    13.43   10/PBC/bptiAndWater_14281
  76.8239    72.28   10/NBC/apoa1AndWater_92224
 179.6466   100.56   10/PBC/apoa1AndWater_92224
 

Script:
protomolLinux examples/water_423/equil298K_01.conf.protomol
namd2Linux examples/water_423/equil298K_01.conf.namd
protomolLinux examples/water_423/equil298K_01.conf.protomol.pbc
namd2Linux examples/water_423/equil298K_01.conf.namd.pbc
protomolLinux examples/bptiAndWater_14281/bpti.conf.protomol
namd2Linux examples/bptiAndWater_14281/bpti.conf.namd
protomolLinux examples/bptiAndWater_14281/bpti.conf.protomol.pbc
namd2Linux examples/bptiAndWater_14281/bpti.conf.namd.pbc
protomolLinux ../TestSuite/apoa1AndWater_92224/apoa1.conf.protomol
namd2Linux ../TestSuite/apoa1AndWater_92224/apoa1.conf.namd
protomolLinux ../TestSuite/apoa1AndWater_92224/apoa1.conf.protomol.pbc
namd2Linux ../TestSuite/apoa1AndWater_92224/apoa1.conf.namd.pbc
 

Did some minor changes in OnAtomPair*.h. I had to add a combined  cutoff force C2/Shift for Coulomb/LJ to use the same switching function as in NAMD2. A first update of MagneticDipole*.[Ch], not verified (Andreas H.'s job)

1)
'restrict' is gone, since it's not part of C/C++ standard. Please feel
free to read more at
http://www.cuj.com/articles/1999/9907/9907d/9907d.htm?topic=articles
I could not see any gain, in some cases it got worse with restrict.

2)
All potential (OneAtomPair*) expect the reciprocal squared
distance. This gave a boost of factor 2 on the IBM using OnAtomPairTwo
(evaluating two potential in the same loop).

3)
Removed several temp variables and const in BSplineMOLLYIntegrator and
HBondMOLLYIntegrator.

4)
The cell manager should run a little faster and it should scale.

5)
The force factory was revisited. It handles properly exact and sloppy
definitions. There was also a tiny problem with lower/uppercase in some
rare cases. Removed an obsolete data member in GenerateIntegrator.

6)
Rewrote PDB, more robust and slim.

7)
Added some statistics for the comparison of forces. It will compute
the average time and the standard deviation.

8)
Some new force for magnetic dipoles were added.

9)
Array.h
http://www.cuj.com/articles/2000/0012/0012c/0012c.htm?topic=articles

Added a more conservative test for the number of atoms in
ExtraTopologyInfo.C. Simulation.C got a more flexible test concerning the
number of atoms in the input files. The flag to abort was set to early in
parseCommandLine.C such that endr in some other function would initiate an
abort. The PSF should be really robust, it synchronizes with !.
ForceFactory parsing problem fixed.

biocomp/TestSuite got some new test case and the old ones where updated.

I did some updates of PME and Ewald such that you can also use a C1 or
Shift switching function. By default it uses Cutoff.

MultiGrid can also be split into two parts: -direct and -smooth

For more information about the input parameters type:
protomol -f

Finally Jan Petter and I were able to fix Nose Hoover NVT Leap-Frog. I fixed also the output of volume, which was undefined when no cell manager was used.

Array.h revisited. resize() will actually ony do re-allocation of memory if the dimensions really differ. This may improve the our code.

MultiGrid, Grid and CubicCellManager using Array<,> do a safe check whether they were able to allocate the memory. In some case when the time step is to big in relation to the dynamics, or when the cell size is to small compared to the simulation box, CubicCellManager will not be able to allocate the requested Array<,>.

SAMPLE got a new column for the volume.

I remove several un-interesting force objects and added some useful two pair forces.

I changed two keywords for the BSDL output such that BSDL related keywords
start with DOBSDL or BSDL. The userguide reflects those changes.

ProtoMol is not able to print out all keywords, default values and
aliases. Just type 'protomol -k' or 'protomol --keywords'.

MultiGrid and CellManager do a safe abort in case they try to request an
array of more memory than the machine can address. This may happen it you
have, say 10 levels for MultiGrid, or when the dimension of the simulation
box is huge and the cell size is tiny.

I brushed up the force object creation and moved it from
GenerateIntegrator into a seperate class: ForceFactory.

Next step will be the IntegratorFactory such that we can add and remove
integrators much easier. The integrator keywords and types will be moved
from GenerateIntegrator into the integrators like it was done with the
forces. The parsing will be much more safe and introducing new integrators
will be pretty easy.

To test all our examples you can go into the exmaples dir and run the
script runAllExamples. It will run all possible configuration for 1 step
and temperature zero (easier to compare, no random velocities).

I removed 48 redundant forces using the algorithm NonbondedSimpleFull,
which are covered by NonbondedCutoff. The binaries shrinked by several
MBytes. Consult your exe using 'protomol -f' to see the supported forces.

I brushed up George's work around for the linking under KCC and did
the same procedure for MIPSpro, adding a small script to link. This
work-around is obsolete for automake 1.5, but this work-around will
generate several warnings when running configure. Now it compiles also on gridur.ntno.no.

stringsutilities's isReal(), isInt(), toReal() and toInt() were rewritten,
based on stdlib (http://www.dinkumware.com/htm_cpl/stdlib.html#strtod)
functions.

Removed dependencies between base add parallel. Parallel will set its own abort function, which is called by Report and others. Fixed the MPI_Win problem for AIX using a wrapper class to avoid <mpi.h> been included inside the header file of Parallel.

Introduced the Paul trap force/potential to simulate icon crystals. Using PaulTrap and Friction as forces we could reproduce some promising results.

The optimized version for SGI under IRIX had some problems with the operator [] in the Array class. Adding a dummy access to the first element fixed the bug. It seems that the IPA does a to 'good' job.

Optimized MultiGrid using references to access arrays for a given level and moved array access into the most outer loop if possible. Gain ca. 2-4x faster. There is still a bug when using production flags on IRIX, but not on Solaris!?!

Could reproduce the error numbers of David H. for C2, C3 and cubic, quintic interpolation, but there seems to be a bug some where in longRangeTerm(). The cell manager for normal boundary conditions was optimized/fixed such it does only consider offsets not bigger than the simulation box. The problem of different error computation made the comparison of David code impossible and the mass for O and H were also different.

Ported successfully the framework to AIX 5.1 on a IBM Power3 machine.

Did the TTP5 homepage, including some new pictures, movies and the 1.05 binaries.

Some first really promising results comparing PME against Ewald:
testsuit/water_3243/water_3243.conf.PBC.PME.Ewald :

PME: alpha=0.422282, V=72324, Rc=7.5, Kc =3.13918, accuracy=1e-06, interpolation=BSpline, order=8.
Ewald: alpha=0.422282, V=72324, Rc=8.80199, Kc (18888)=3.13918, accuracy=1e-06, misses=3.58429%.

error : F2 5.14152e-05, Fmax 6.66113e-05, PE 2.43768e-05
error : F2 5.12943e-05, Fmax 6.75503e-05, PE 2.42813e-05
error : F2 5.13815e-05, Fmax 6.73303e-05, PE 2.41761e-05
error : F2 5.13703e-05, Fmax 6.81136e-05, PE 2.3929e-05
error : F2 5.13466e-05, Fmax 6.53534e-05, PE 2.4129e-05
error : F2 5.13122e-05, Fmax 6.57097e-05, PE 2.40598e-05
error : F2 5.14361e-05, Fmax 7.52524e-05, PE 2.38136e-05
error : F2 5.13735e-05, Fmax 7.14231e-05, PE 2.38519e-05
error : F2 5.14119e-05, Fmax 7.30841e-05, PE 2.37703e-05
error : F2 5.13692e-05, Fmax 7.43994e-05, PE 2.39662e-05
error : F2 5.14211e-05, Fmax 7.59528e-05, PE 2.43565e-05
Timing 'PMEwald' (0) [s]: user time 14.39
Timing 'FullEwald' (0) [s]: user time 194.6

1) I added a new application namd2protomol, but just the interface.
 

2) I changed the application directories from app-* to *-app.
 

3) The Lagrange interpolation was originally a Hermite interpolation.
Lagrange has for the time being a empty body, where Hermite is the old
code from Lagrange. Additionally the quintic order *Hermite) was added.
 

4) The force definitions for PME, FullEwald and MultiGrid changed such
they are consistent with all other non-bonded forces:

  Coulomb -algorithm FullEwald -real -reciprocal -correction
    -alpha <real=-1>
    -accuracy <real=1e-05>

  Coulomb -algorithm MultiGrid -interpolation BSpline -kernel C1
    -toplevelgrid <int> <int> <int>
    -levels <int>
    -s <real>
    -order <int=4>
    -ratio <int=2>

  Coulomb -algorithm PMEwald -real -correction -interpolation BSpline
    -gridsize <int> <int> <int>
    -cutoff <real>
    -order <int=4>
    -accuracy <real=1e-06>
    -alpha <real=-1>
 

Presented a status report of ProtoMol at the TTP5 meeting in Tromsoe.

ProtoMol.html: Added two new movies.
 

CompareForce: ProtoMol is now able to compare two forces and energies on the fly. It
computes the relative error. To compare forces just add the keyword
compare after the keyword force. Ex. Ewald vs Direct:

force compare FullEwald
force compare Coulomb -algorithm NonbondedSimpleFull

NonbondedSimpleFull is the reference force, where FullEwald is the
approximation.(http://www.cs.ut.ee/~toomas_l/linalg/lin1/node9.html Def
1.5.5).

Additionally it prints also the timing for each force. You may also
compare more than one pair of forces in the same run , Ex. Coulomb and LJ.
Only the approximated force contributions are added to the system.

CompareForce concept is based on polymorph inheritance.
 

EnergyStructure: I moved EnergyStructure to framework/topology and added two sum
methods:

sum(topo,positions)
sumPE()
 

NonbondedFullEwaldSystemForce & NonbondedPMEwaldSystemForce: The argument and keyword -j <real> are removed for PBC since it does not
make sense.
 

PROTOMOL is an object-oriented component based framework for molecular
dynamics simulations.  The framework  supports the CHARMM 19 and 28a2 force
fields and is able to process PDB, PSF,  XYZ  and  DCD trajectory files. It
is designed for high flexibility, easy extendibility and maintenance, and
high performance demands, including parallelization.  The technique of
multiple time-stepping has been used to improve long-term efficiency, and
the use of fast electrostatic force evaluation algorithms like plain Ewald ,
Particle Mesh Ewald , and Multigrid summation   further improves
performance.   Longer time steps are possible using  MOLLY,  Langevin Molly
and Hybrid Monte Carlo ,  Nose-Hoover, and Langevin integrators. In
addition, PROTOMOL has been designed to interact with VMD, a visualization
engine developed by the University of Illinois that is used for displaying
large biomolecular systems in three dimensions.  PROTOMOL is free
distributed software, and the source code is included.

http://www.ks.uiuc.edu/Development/biosoftdb/biosoft.cgi?&category=2

PME and Ewald can be used also for the normal boundary conditions. It will
just add space around the actual simulation box, which is specified by '-j
<real>'. The default value is 3, which means you get one shell around the
simulation box with thickness of the actual simulation box. 5 would
add 2 shells.

Multigrid should work.

The cell manager works now in O(N) and will abort in case the system
expands to infinit (when timestep is to big).

Marc is working on BBKImpulseIntegrator and LangevinImpulseIntegrator.

I added a new application analyzer, which reads frocefiles and compares
them and computes the relative error. It's also an example for future
applications. Be aware that there is already a cmd analyzer ;-)

r : 1 [AAm]    = 1e-10 [m]
t : 1 [s']     = 1e+15 * sqrt((1.660ee-27 * 1e-20 * 6.022035e+23)/4184)[fs]
v : 1 [AAm/s'] = 1e+15 * sqrt((1.660ee-27 * 1e-20 * 6.022035e+23)/4184)[m/fs]
F : 1 [kcal/mol AAm] = 4184 / (6.022045e+23 * 1e-10)[N]
E : 1 [kcal/mol]     = 4184 / 6.022045e+23 [J]
m : 1 [amu]          = 1.6605e-27 [kg]
e : 1 [e]            = 1.602e-19 [C]

TTP WorkShop last week:

1) Drop-drop and drop-homophase coalescence
We were able to convert their input format to PSF/PAR/PDB/XYZ and run a
simulation with those droplets:
http://www.chembio.ntnu.no/~bhaf/forskning/forskning.html#coalescence

Here I see some new requiremnts ... or more a HowToDo document to explain
how you add code at the right place to generate user specific ouput.

2) Magnetic Dipole Force

I found only this (but will get some info ...):

In 1979, Professor John Ugelstad of Norway achieved something under normal
gravity otherwise only achieved by NASA in the weightless conditions of
space - the making of uniform polystyrene spheres of exactly the same
size. Prof. Ugelstad and his colleagues made these spheres magnetisable
and, moreover, superparamagnetic meaning they are only magnetic in a
magnetic field. Due to this property, the particles can easily be
resuspended when the magnetic field is removed. This discovery has
revolutionised the isolation and separation of many biological materials.
For instance, the attachment of target-specific antibodies to the surface
of the beads allows capture and isolation of intact cells directly from a
complex suspension, such as blood. This is all
accomplished under the influence of a simple magnetic field without the
need of column techniques or centrifugation.

Friday     2001.10.19

Major update and adding new functionalities contributed by Andreas Hellesoey:

- PME works on SGI and other systems using FFT. ProtoMol matchs NAMD2.2.
  Interpolation order and cutoff can be specified.
- gccreturn.h is history.
- Removed some redundant force objects (NonbondedFullSystemForce<> using
  NormalBoundaryCondtions) makes the exe smaller.
- Removed templates for bonded forces since we do not use them,
  makes exe smaller (removed also the bonded base class files).
- The compilation is some sort faster.
- '-DTEMPLATE_IN_HEADER -DEXPORT=' flags are history.
- Added a substitute for the SGI FFT -lcomplib.sgimath, a translation of
  Fortran code (http://fy.chalmers.se/~tfylb/ZFFT.html). It's small &
  fast, but does not support all sort of transformation length. NB: It
  creates 1K warnings ...
- NoseNVT fixed.
- Updated the doc (NoseNVT and PMEwald)
- Update the ProtoMol homepage and added ProtoMol Alpha 1.03 binaries.
- Full optimization for IRIX MPI works again.
- Added new class for the B-splines used in PME (NAMD like, but nicer).
- Added a new base class for ExtendedFroce and SystemForce to make them
  transparent to the user and simplify the parsing.
- MagneticDipoleForce is a new potential, which is non-radial.
- The interfaces for doOneAtom was extended by the diff vector for
  non-radial forces.
- Implemented a first example of a 1-body extended force depending on the
  velocity (FrictionExtendedForce).
- If possible I made the parsing methods static.
- 'protomol -f' prints out all supported forces with their keywords,
  arguments, type and default values.
- protomol changes to the directory where the config file/prefix file is
  located.
- protomol writes back all integrators and forces, from their internal
  settings.
- Non memory leaks, beside from the SGI-FFT-lib.
- Corrected 'firststep' of the restart config files.

- Major rework of GenerateIntegrator + small bug fixes (DCDWriter, reading
  XYZ pos and no PDB)
- All examples checked

Each possible/ desired force type is instantiated and registered
(framework/frontend/registerAllExemplars*.C) only once without any
arguments/ parameters. When ever a force is needed a copy is made and
its arguments/ parameters are initialized (make(vector<VarVal>)). Each force
provides a method (getParameters()), which returns a vector of
keyword, type, value and default value. The vector passed to make() must
correspond with the VarVal of the to vector of getParameters(). See
SystemForce.h

I did also rework the parsing of forces such that it accepts a sloppy
definition and an exact one. The sloppy version is as it was before, where
the order and the exact number of parameters (e.g same cutoff for
algorithm and the switching function, or same switching function for LJ &
Coulomb) do not matter. The exact version is useful e.g. HBondMolly,
where values with the same keyword may differ (e.g. C1 and C2 switching
functions). The parsing try first to find an exact match to resolve the
correct force type with its templates. It requires that the input for the
force type matches the string (getIdString()) of one registered force
object. If more keywords match a force type using sloppy parsing than
under exact parsing, then sloppy definition is preferred. For the
arguments/ parameters yields the same. An exact parsing of arguments/
parameters presumes an exact fore type definition.

With this new approach we also get the write back of forces for free,
and the parsing gets more user friendly when an error occurs. All
possible forces with their arguments are/ can be printed out.

To make the compilation more comfortable I spread the instantiation of
force objects over several files. You will see more progress and the
memory required to compile such be less. When developing new forces
you may also comment out some force objects to speed up the compilation.
E.g. comment out '#define USE_NONBONDEDFULLSYSTEMFORCE' in
framework/frontend/registerAllExemplarsFull.C.

I fixed all examples according to the new keywords, but there are
still two examples not working due to a reader (chram28a) problem
(2mlt.conf, enkphalin_75). I fixed also the water_6 example, where I
had to rework the part of ExtraTopologyInfo when no PDB file exists,
only a XYZ one. Now it is possible to run it only with psf & xyz.

It seems that g++ does not like some special order of different
types under production compilation. By changing the order of
initialize the bus error disappears. I comment out the file option
'std::ios::app' suppressed any output file.

To make the parsing not case sensitive we could use our own string
comperator (ltstrNocase in stringutilities.h) for all maps using
string as key, but the compilation time for production code on SUN &
g++ would increase with ~10'.

Ran the new front end with periodic boundary conditions. No problem!

Fixed the segv for the production code on SUN for g++. I changed the order of arguments of the integrator constructors such that we start to pass atomic types, objects by value and last all pointers. It seems to help g++ ...

Removed virtual warmUP(){} in StandardIntegrator.

Added a new method equalNocase(const string&,const string) for no-case sensitive compare. I replaced all equalBeginNocase in GenerateIntagrator to avoid some really strange parsing behavior, e.g. '-algorithm Nonbond' would match using equalBeginNocase.

Yes! A perfect match of energies with the old front end for alanin up to 7000 steps! Note that you have to set the temperatur to 0 or make sure that you generate the same intial velocities.

Finally found the problem casuing the wrong nonbonded energies:

The setting of the exclusionType was not done correctly. Use 'VarVal(ExclusionType(ExclusionType::ONE3))',
ExclusionType::ONE3 is implicitly converted to int ...
 

When comparing, please check that you really have the same configuration(integrators, forces, switchingfunctions, temperature, etc. ...). I fixed the alanin example: SWC1 -> FSWC1.
 

stringutilities:
----------------
equalEnd(...)
equalEndNocase(...)
equalBegin(...)
equalBeginNocase(...)

I replaced all cmp_nocase(string,string,int) by equalBeginNocase(...),
except in Par.C, since it seems that it does not work.
Better to let the function figure out the number of chars to compare.
equalBegin*() compares the first min(s1.size(),s2.size()) chars, where
equalEnd*() compares the last min(s1.size(),s2.size()) chars. As I can see
there is no need to use cmp_nocase with number.
 

ExclusionType:
--------------
Using AbstractEnumType as base class to write advanced enum types you get
some very handy functionality. As you can see below, you do not need to
know what ExclusionType takes as values. If you want to add a new
exclusion type, just add it in ExclusionType.[Ch] ... the whole idea is
based on conversion operators. There is only one but (as I know): it's
case sensitive ... assuming that the string values of ExclusionType are
lower case just pass a lowercase string.

--------------------------------------------------------------------------
else if (myDataTypes[index] == "exclusion_type"){
  string temp;
  configFileStream >> temp;
  if (cmp_nocase(temp, "none", 4))
    (myData->lookup(index)).setExclusion_Type(ExclusionType::NONE);
  else if (cmp_nocase(temp, "1-2", 3))
    (myData->lookup(index)).setExclusion_Type(ExclusionType::ONE2);
  else if (cmp_nocase(temp, "1-3", 3))
    (myData->lookup(index)).setExclusion_Type(ExclusionType::ONE3);
  else if (cmp_nocase(temp, "1-4", 3))
    (myData->lookup(index)).setExclusion_Type(ExclusionType::ONE4);
  else if (cmp_nocase(temp, "scaled1-4", 9))
    (myData->lookup(index)).setExclusion_Type(ExclusionType::ONE4MODIFIED);
  else{
      report << error << "[SimulationConfiguration::configRead]"
             << temp << " is not a valid exclusion type in"
             << " the configuration file.  "
             << "Valid types are: none, 1-2, 1-3, 1-4, or"
             << " scaled1-4.  Exiting." << endr;
  }
}

--------------------------------------------------------------------------
else if (myDataTypes[index] == "exclusion_type"){
        string temp;
        configFileStream >> temp;
        ExclusionType exlucsionType(lowercase(temp));
        if(exlucsionType != ExclusionType::UNDEFINED){
          (myData->lookup(index)).setExclusion_Type(exlucsionType);
        }
        else {
          report << error << "[SimulationConfiguration::configRead]"
                 << temp << " is not a valid exclusion type in"
                 << " the configuration file.  "
                 << "Valid types are: "<<ExclusionType::getPossibleValues(", ")
                 << ".  Exiting." << endr;
        }
}
--------------------------------------------------------------------------

myDataTypes:
------------
I did some tests to replace the type of SimulationConfiguration::myDataTypes
by map<string, VarValType>, but then it does not work because of the type
"integrator". To extend VarVal with INTEGRATOR is not really what we want
(?). What I did so far, was to replace all "string", "int", etc. by the
strings of VarValType. I'm not that happy how the types are assigned to
myDataTypes, but I can not figure out how I've to overload or implement
the converions operators ... the advantage is that the correct/valid
strings are assigned, e.g. avoiding 'strign' or like that. If you do not
like/ not agree you can simply change it back with replace.
 

VarVal:
-------
I removed VarVal::writeData since we can get a string of VarVal by
VarVal::getString() or VarVal::toString(). We should also allow some more
conversion between types, especially  any type <-> string. We would not
need to to the explicit conversions ... 'let papa do the rest'.

SimulationConfiguration::myGlobals:
-----------------------------------
I removed SimulationConfiguration::myGlobals. It retrieves the boundary
conditions and the cell manager at the right position and place
(Simulation). It's only local stuff and we still have the (string-)values
in myData if later needed.

pdbWriter():
------------
The pdbWriter is not 100% correct ... the coords are wrong. I fixed the
writer such that it writes 80 chars. You may fix the rest. My suggestion
is to use temp string for one line and use toString rather use
stringstreams directly. As you can see you can write down one column in
more or less one line. One minor note: I commented out the endings in
pdbWriter and psfWriter since they are already added before the call.

Now I get a segfault after the first timestep ...
 

Debugging ...zzzZZZzzz... more fixing code to be more elegant ...

I figured out why Simulation crashes. Your definition of member data
in PSF and Par and the use of integrators does not work together. Every
time you call getAtoms(), getBonds(), etc. you will get a new copy ... the
for loop will not work properly since .begin() and .end() are not of the
same object.

I suggest to make all data members access by get*Pointer()->get*() public
since I do not really see reason to make a copy every time to what to
access it. I may use pointers to define the data members, but you will
still pass a pointer such data the data can be modified.

I attach my changes so far. You may take a look at XYZ. I did some
changes, but I think you should take a close look ... it should look like
PSF, Par, etc.. Is it appropriated to use pointers for myXYZBlock? I do
not see any new nor do I find any destructor.

I did a minor fix to cmp_nocase().

To use pointer for Coordinates inside PDBElement seems to me not really
appropriated. The constructor use member initiation list such you can
define a const object.

I guess there was an error cause to wrong indexing. The atoms start with
0, where they start with 1 in the input files ... I reworked atoms, bonds,
dihedrals and impropers.

Did some small performance test of Array<>, and compared with f90 on the Origin2000. Compiling the test case without assert's, Array<> ca. 2-3 slower than the FORTRAN executable. Most important is who elements are access in nested loops. Wrong nesting can slow down up to a factor 30-40! This yields both for C++ and f90.

Started to add new system force using multigrid methods. Implemented a Lagrange class for the inter-/anterpolation. We may re-think about how such interpolation methods should be inserted into ProtoMol

Helping Tom and Trevor to debug ...zzzZZZzzz...
 

The energies of ProtoMol (EquilibriumMOLLY 4cyclelength + Leapfrog 1.0fs) match those from NAMD2 for equil298K_01, but we can see some differences after 1000 [fs].

The walltime for ProtoMol running EquilibriumMOLLY/Leapfrog 2000 [fs] is 90 [s] vs. NAMD2 54 [s]. I guess there is still improvement possible.

The energies of ProtoMol (Leapfrog 1.0fs) match those from NAMD2 for equil298K_01, but we can see some differences after 1000 [fs].

The walltime for ProtoMol running Leapfrog 2000 [fs] is 62 [s] vs. NAMD2 54 [s].

Memory usage (top). Good to see that we allocate what we really need/use:
     PID        PGRP USERNAME PRI  SIZE   RES STATE    TIME WCPU% CPU% COMMAND
   6374400    6374400 matthey   20   33M 6368K run/97   0:44 81.4 99.67 namd2
   6407632    6407632 matthey   20 6816K 4432K ru/117   0:43 90.1 99.70 protomo
 

ProtoMol:

 TS   Angle      Bond          Coulomb  VdW              Kin       Total      Temp

   0 121.42818   76.21106   -1805.24376 226.01204     0.00000   -1381.59248   0.00000 
 100  82.37828   69.43713   -2180.04433 377.66611   262.44232   -1388.12048 208.14339 
 200  99.63987   70.82989   -2103.59479 308.80243   238.16935   -1386.15326 188.89246 
 300  98.80270   77.35538   -2135.38167 308.17179   264.85648   -1386.19531 210.05807 
 400  92.93471   78.17219   -2156.85586 331.92604   266.72687   -1387.09605 211.54147 
 500  95.77862   74.46579   -2146.98236 319.04947   270.97403   -1386.71445 214.90991 
 600  98.13376   83.75192   -2169.61434 333.14127   268.01478   -1386.57262 212.56292 
 700  97.09825   75.66248   -2110.79272 286.24480   264.99492   -1386.79227 210.16786 
 800  96.77886   74.04126   -2120.85556 301.68168   261.56765   -1386.78611 207.44969 
 900 102.74433   84.26903   -2159.05949 329.36597   255.95236   -1386.72780 202.99619 
1000 101.32859   80.19226   -2136.17258 303.42246   264.52723   -1386.70205 209.79694 
1100 103.61514   67.12995   -2139.76204 310.96207   271.14706   -1386.90781 215.04713 
1200 104.17008   81.59456   -2148.97060 299.29573   277.02385   -1386.88638 219.70802 
1300 114.37683  105.99985   -2166.61733 311.94094   247.86870   -1386.43100 196.58503 
1400 107.61732   81.59097   -2178.32416 328.97391   272.56814   -1387.57383 216.17419 
1500 115.09214   73.25639   -2118.20763 290.55711   252.87076   -1386.43123 200.55218 
1600 112.73964   93.27955   -2171.26582 307.76209   270.31793   -1387.16661 214.38955 
1700 119.48309   93.83196   -2174.92292 297.14384   277.78843   -1386.67558 220.31442 
1800 117.23035   70.62509   -2158.87230 307.28310   276.84663   -1386.88714 219.56747 
1900 117.17593   78.14015   -2146.83178 286.25681   278.27652   -1386.98238 220.70152 
2000 119.91404   82.33398   -2164.76961 300.03900   275.40045   -1387.08214 218.42051

NAMD2:
 TS   Angle      Bond          Coulomb  VdW              Kin       Total      Temp

   0 121.4282    76.2111    -1805.2437  226.0120      0.0000    -1381.5924    0.0000    
 100  82.3782    69.4372    -2180.0442  377.6661    262.4423    -1388.1204  208.6366  
 200  99.6397    70.8299    -2103.5947  308.8024    238.1695    -1386.1532  189.3402  
 300  98.8027    77.3553    -2135.3815  308.1718    264.8565    -1386.1953  210.5558  
 400  92.9346    78.1721    -2156.8557  331.9260    266.7269    -1387.0960  212.0428  
 500  95.7785    74.4658    -2146.9822  319.0494    270.9741    -1386.7144  215.4192  
 600  98.1334    83.7519    -2169.6141  333.1412    268.0149    -1386.5726  213.0668  
 700  97.0983    75.6625    -2110.7928  286.2444    264.9953    -1386.7922  210.6662  
 800  96.7787    74.0413    -2120.8599  301.6845    261.5693    -1386.7861  207.9426  
 900 102.7439    84.2665    -2159.0408  329.3549    255.9480    -1386.7275  203.4737  
1000 101.3272    80.2011    -2136.1900  303.4160    264.5429    -1386.7029  210.3065  
1100 103.6146    67.1527    -2139.7295  310.8176    271.2344    -1386.9102  215.6262  
1200 104.1137    81.5707    -2148.4123  299.1808    276.6624    -1386.8848  219.9413  
1300 114.5995   105.8455    -2169.6823  311.8970    250.8992    -1386.4411  199.4601  
1400 108.2804    79.8893    -2176.8184  328.5486    272.5437    -1387.5563  216.6670  
1500 113.2433    70.4257    -2126.6604  294.2870    262.2653    -1386.4391  208.4959  
1600 108.4011    91.8701    -2152.4721  291.7157    273.2904    -1387.1948  217.2606  
1700 115.1490    92.6138    -2179.6727  312.1955    272.9166    -1386.7977  216.9635  
1800 115.9492    85.1862    -2159.2748  295.9645    275.4964    -1386.6785  219.0143  
1900 108.8865    69.1596    -2143.4247  300.6773    277.4422    -1387.2591  220.5612  
2000 111.4124    92.7068    -2161.7919  303.8573    266.8553    -1386.9600  212.1449

ProtoMol:

initial_t 0.0    
nsteps 2000
cell_size 6.5
outputf 20
init_temp 0

boundaryConditions Normal
cellManager Cubic

Integrator {


  level 0 Leapfrog {
        timestep 1.0
    force Bond, Angle 
    force Coulomb 
          -algorithm NonbondedCutoff
          -switchingFunction Shift
          -cutoff 6.5
          -switchon 0.1
    force LennardJones
          -algorithm NonbondedCutoff    
          -switchingFunction FSWC1
          -cutoff 6.5
          -switchon 0.1
  }

}



NAMD2:

structure       equil298K_01.psf

parameters      equil298K_01.par
exclude         1-3
1-4scaling      1.0

switching       on
switchdist      0.1
cutoff          6.5
margin          0
stepspercycle   4

molly off

timestep 1.0

outputenergies  20
outputtiming    20
binaryoutput    no

coordinates     equil298K_01.pdb

outputname      equil298K_01.namd.output

temperature     0
numsteps        2000

Total energy of PMEwald still not matching FullEwald.

ProtoMol matches (+1%) the performance of NAMD2 (own compilation) for equil298K_01.

TTP meeting. Presentation of ProtoMol.

switchingDeriv in ShiftSwitchingFunction was wrongly computed, factor sqrt(distSquared) removed. Using ShiftSwitchingFuction we get a perfect match of energies with NAMD2 numbers, over several time steps! Could not reproduce any test run case under NAMD2 using equilibrium.

Moved setForceGroup(); back/up to MTSIntegrator as it was previous. Changed back GenerateIntegrator such that forceGroup never points to a zero pointer for the instantiation of the integrators, but it may be empty.

Reworked HMCIntegrator::initialize(...) such that it does execute the warm up. Removed warmUp(), also the hack in protomol.C calling warmUp().

Ran Virgils PME code successfully. Remains only to optimize and verify it ...

ReadVelocFile(...) produces a memory leak: line 'double* s  = new double[numberAtoms]();'

I added ShiftSwitchingFunction as implemented in NAMD2 and could reproduce the same Coulomb energy number for the first time step only. ShiftSwitchingFunction maybe wrong?

Verifying EquilibriumMOLLY against NAMD2. ProtoMol does produce the same energies only using NonbondedSimpleFull and Leapfrog. When using cutoff the Coulomb energy does not match NAMD2 numbers. No explanation for the moment. Next step is to download the NAMD2 src and do one to one force compare.

NAMD 2.2 compiled and linked. Number of interactions is the same for NAMD and ProtoMol using cutoff and Leapfrog, but the Coulomb energy is different for each interaction ... some factor.

New data member for AtomType.

AtomType gets a new data member symboleName, which holds the atom unity name. E.g. needed to figure out all H atoms regardless their different representation. Assigned to Trevor.
 

Tryied to figure out the compiling problem on Solaris 5.8 with KCC.

Tom:

> I noticed while working on Stones that GenerateIntegrator does not
> compile with KCC on Solaris 8.  The error seems pretty bizare:
>
> KCC -DHAVE_CONFIG_H -I. -I. -I../.. -I../../framework/autogenerate
> -I../../framework/base -I../../framework/topology
> -I../../framework/integrators -I../../framework/parallel
> -I../../framework/forces -I../../framework/frontend -I../../framework/imd
> +K0   --one_per   --restrict -DVECTORTYPE_NVDOUBLE -DTEMPLATE_IN_HEADER
> -DEXPORT= -DMPI_NO_CPPBIND -c GenerateIntegrator.C
> "/usr/include/unistd.h", line 382: error: function
>           "rename(const char *, const char *) C" conflicts with
>           using-declaration of function
>           "std::rename(const char *, const char *) C"
>   extern int rename(const char *, const char *);

I think this may be one good reason why we should not use 'using namespace std;'. I did a grep over all includes of GenerateIntegrator.C on stones:

/afs/nd.edu/sun4x_57/usr/local/src/kai/kai34g/KCC_BASE/include/fstream:1:    using namespace std;
/afs/nd.edu/sun4x_57/usr/local/src/kai/kai34g/KCC_BASE/include/locale:93:using namespace std;
/afs/nd.edu/sun4x_57/usr/local/src/kai/kai34g/KCC_BASE/include/mcompile.h:7:#define MSIPL_USING_STD using
namespace std;
 

But I could not reproduce the error on IRIX ... and the 'using namespace
std;' are local, inside the namespace __kcc. What I can say is that stones
and guadalupe include the same files under KCC_BASE and
/usr/include/unistd.h should define the same rename for 5.7
and 5.8.

Tried to compare the preprocessed files (-E option) of
GenerateIntegrator.C for 5.8 and 5.7  ...

GI57.C:extern int       rename(const char *, const char *);
GI57.C:using ::rename;
GI57.C:extern int rename(const char *, const char *);

GI58.C:extern int       rename(const char *, const char *);
GI58.C:using std::rename;
GI58.C:using ::rename;
GI58.C:extern int rename(const char *, const char *);

I guess 'using std::rename;' does cause the problem, inside
/usr/include/unistd.h.