Fall 2007
*********************************************************************************************************
DOWNLOADS:
*********************************************************************************************************
*********************************************************************************************************
PAPERS:
ProtoMol How-To (PDF)
Current Thesis (PDF)
*********************************************************************************************************
*********************************************************************************************************
NOTES:
TAU - PERFORMANCE PROFILER THAT WORKS WITH PYTHON
Paul's Alanine Trajectories are here.
*NEW* TARGET CONFERENCE FOR FTSM here.
CHECK OUT BIOPERL: http://www.bioperl.org/wiki/Main_Page
PDYNAMO
SUMMARIZING MY WORK AS 3 PSEs
DR. IZAGUIRRE NUMBER 650-724-9826
*********************************************************************************************************
*********************************************************************************************************
USEFUL PERSONAL LINKS:
ND THESIS: STYLEGUIDELINES
REFEREEING PAPERS
SPARC AUTHOR RIGHTS
LCLS NORMAL MODE SITE
REMOVE .svn FILES RECURSIVELY: find . -name .svn -exec /bin/rm -r {} \;
GET DETAILED ARCHITECTURE INFO: cat /proc/cpuinfo | more
*********************************************************************************************************
*********************************************************************************************************
TODO:
*********************************************************************************************************
Sunday, Dec 23
I think at this point I have everything I want syntactically from MDL written for the JCP paper.
The next step will be to talk with Chris and better understand the following methods:
1. NosePoincare
2. RMT
3. Takahashi
4.Parallel In Time
I'll spend the rest of the day on FTSM.
Friday, Dec 21
Peter wrote me this morning and appeared happy after the bug fixes from yesterday.
One thing I'll have to make sure of now is that scpenergy-app is also accurate (I only changed
tools/scpenergy).
Can't print anything, but I do need to clean up the Parallel In Time simulations in MDL, they
're still using old notation..
I uncovered a huge bug with FTSM this afternoon! The target (phi, psi) values were not set properly
for the constrained dynamics runs. They were being set to the replicas one above the current.
The string no longer shifts upward after this change.
Results - the solvated simulation works a little bit better, but smoothing (like always) screws that one
up. This is an attempt to run the same system that EVH and Luca ran, and they started their string
with eight points between (-40, 130) and (-40, -45). The fact that all eight points are almost in a straight
vertical line, somehow the smoothing routine does not like it. If I turn off smoothing, results look better
initially, albeit the string of course eventually evolves to a porcupine.
The unsolvated is still all over the map; however I just have no idea which parameters are accurate for
that type of system.
Wednesday, Dec 19
Finally finished with the grading!
Now looking at the scp woes again...
Okay, I am narrowing down the differences in CoulombSCPISMForce. It seems like all
of the SCPISM-related parameters are different.
Tuesday, Dec 18
Wrote a lot more with the JCP this morning.
I'm becoming more and more convinced that I'm really going to need to merge the linear code
with CVS protomol. Keeping linprotomol up-to-date with the updates to protomol is just too much
overhead.
And eventually, it may even be convenient to have mdl/ within the protomol repository.
You have to install ProtoMol anyway to use mdl, because you need to be able to compile the shared
objects locally.
There are going to be two issues here:
(1) I can only overload operator[] for Vector3DBlock to return a value but not a reference.
Since Vector3DBlock will now wrap a C array, if you index I can return a newly created Vector3D object
but I cannot reference data in the C array. So we'd have this situation:
Vector3DBlock v(24);
Vector3D v3d = v[12]; // Legal
v[14] = v3d; // Illegal
(2) Some of the code from protomol will be incompatible with SWIG, requiring the use of #defines.
Anyway, looking again at Peter's troubles..
Monday, Dec 17
Fixed a couple small errors this morning:
1. Made quartic switching, parameters take effect if doscpism = 3.
2. Fixed the cutoff in scpenergy-app and tools/scpenergy to be a Real (not an int).
I also believe I fixed Pauls' coinfigure issues on helios with icc, but I'm going to try it everywhere just to make sure.
Tests I should run (hiliting as they work):
Machine Compiler MPICH LAPACK
----------- ----------- --------- -----------
helios icpc N N
helios icpc N Y
helios gcc N N
helios gcc N Y
iss gcc N N
iss gcc N Y
iss gcc Y N
iss gcc Y Y
opteron icpc N N
opteron icpc N Y
opteron icpc Y N
opteron icpc Y Y
opteron gcc N N
opteron gcc N Y
opteron gcc Y N
opteron gcc Y Y
Okay, the two gcc mpich tests failed on opteron (it means they'd probably fail anywhere...)
Wait, no I'm using an incompatible simtklapack it appears.
But wait, how did the other ones work? No it can't be.
Okay gcc-mpich works. It's just gcc-mpich with lapack now that screws up.
Fixed, 32-bit mpicxx does not work on hpcc. 64-bit does, but for that
I needed 64-bit Lapack. Okay, we're good.
Now writing...
Sunday, Dec 16
Going to try to nail these new parameters this morning, and to instantiate our CoulombBornRadiiForce
behind the scenes...
Here's the output I have to match:
Step : 0, Time : 0.000 [fs], TE : -292.0603 [kcal/mol], T : 462.4952 [K], V : 202.02 [AA^3]
Step : 1, Time : 1.000 [fs], TE : -296.9108 [kcal/mol], T : 344.5921 [K], V : 202.13 [AA^3]
Step : 2, Time : 2.000 [fs], TE : -296.8111 [kcal/mol], T : 307.1401 [K], V : 201.29 [AA^3]
Step : 3, Time : 3.000 [fs], TE : -295.4440 [kcal/mol], T : 302.4017 [K], V : 199.63 [AA^3]
Step : 4, Time : 4.000 [fs], TE : -295.9702 [kcal/mol], T : 270.1764 [K], V : 198.33 [AA^3]
Step : 5, Time : 5.000 [fs], TE : -298.6648 [kcal/mol], T : 220.4174 [K], V : 197.71 [AA^3]
Step : 6, Time : 6.000 [fs], TE : -298.7586 [kcal/mol], T : 236.7842 [K], V : 198.81 [AA^3]
Step : 7, Time : 7.000 [fs], TE : -298.4081 [kcal/mol], T : 289.2521 [K], V : 200.16 [AA^3]
Step : 8, Time : 8.000 [fs], TE : -301.5930 [kcal/mol], T : 297.0191 [K], V : 202.66 [AA^3]
Step : 9, Time : 9.000 [fs], TE : -304.0111 [kcal/mol], T : 293.6348 [K], V : 204.54 [AA^3]
Step : 10, Time : 10.000 [fs], TE : -304.4633 [kcal/mol], T : 290.0030 [K], V : 206.15 [AA^3]
Let's see...I will have to be careful that if the user DOES specify CoulombBornRadii, that I don't insert the same force twice.
YAY! It works!!!
Friday, Dec 14
Thierry's e-mail convinced me of the merits of switching functions for the Full algorithm.
I do think those should stay in MDL. However, I don't think the SimpleFull ones should be there
especially if their sole purpose in ProtoMol was to validate that Cutoff worked with
the same switching functions. Those have bsince been tested, and in reality no one would apply
a SimpleFull algorithm with a switching function. I'm going to remove them.
Okay, now onward to Hassan's parameters and the SCPQuartic.
I have put the new parameters into ProtoMol - only two changes from Hassan's new ones; the hbond_factor
for H and HC atomtypes. This is such a small change that I did not build a second table.
Somehow it did not affect the results even though we have
a couple of H atomtypes in the simulation.
Thursday, Dec 13
Cleaned up the factory a bit more; also am improving the structure of the Python prototyped forces, giving them their own
directory, etc.
Now the one I am committing (I forgot about this) uses the libraries which approximate derivatives, and so
the force calculation is a bit off.
I think I'm going to change it, to commit one that is accurate.
I also have not yet prototyped Python forces which apply the correct boundary conditions. That needs to be done as well.
Actually I'm going to commit both - since this is the best 'smart' implementation I have thus far.
Actually, in a weird way the "dumb"version actually turned out to be simpler!
Biggest reason: we needed to set up the dihedral function of the positions in Python,
sot that the derivative computation could view the symbol 'x' and execute the chain rule.
Writing more of the JCP paper.
I'm just realizing something... why do we have switching functions for SimpleFull and Full?
They don't make sense if we're not using a cutoff. Right??
Wednesday, Dec 12
The force factory works great now. All of the method handles for the SWIG-wrapped object
constructors are registered in a Python dictionary in the ForceFactory constructors, and all those HUGE
'if' statements are now gone. I've shrunk the file by about 25%.
The remaining problem is this MultiGrid business.
Clearly some parameters are not being set properly. I'm doing a run with ProtoMol using alanalogous
parameter values. Here is the difference:
CORRECT: MultiGrid: s=10, ratio=2, h=2.45481,2.11413,2.29575, Hermite 6th order C3, levels=(16,16,16),(8,8,8), finest grid=(-12.491,-11.799,-15.251)-(26.786,22.027,21.481), Gc(9,10,9).
MINE: MultiGrid: s=10, ratio=2, h=3,3,3, Hermite 6th order C3, levels=(21,19,19),(15,14,14), finest grid=(-21,-21,-24)-(39,33,30), Gc(7,7,7).
Looking into this...
Fixed it. H should not be manually set for PBC, but should use default values.
Oh yeah! Functional MultiGrid!
This factory is looking excellent. I reduced code duplication by handling parameter setting in a single function
dependent upon the algorithm and switching function(s). Now we're at about 600 lines (used to be 1200).
Tuesday, Dec 11
This force factory is *really* looking a lot nicer through the use of Python dictionaries.
All these huge 'if' clauses are going to be gone.
The only thing I'm having a bit of an issue with are the extra parameters... sometimes they should be used
and sometimes they should not. I'm thinking of adding a second dictionary for these.
Another difficulty is that depending on the type of algorithm for nonbonded calculations, the number of templated
variables are going to vary.
This is looking better - seems that MultiGrid was broken though with more recent updates from PM. :(
I'll have to debug this again.
Monday, Dec 10
FINALLY got my AFS back, made the changes to scpenergy-app.
It was just as well, my office hours were flooded today and things are busy with the course now at the
end of the semester.
Fixing up the force factory...
Thursday, Dec 6
Wrapped up the F@H units early this afternoon.
Also fixed another bug in the configure.in script. Hopefully, the last.
Now, I need to produce slides etc. for class tomorrow, some new special applications of
recursive fractals that I need to become familiar with.
Plus, need to make a studyguide for the final.
Wednesday, Dec 5
More on the units today... finished NormalMode test.
Tuesday, Dec 4
Working on the work units for Joseph, I want to cover these three examples:
1. STS
2. Normal Modes
3. SCPISM
Normal modes will cover the MTS case with NormModeInt and NormModeMin.
I'll need to come up with a full set of parameters to vary.
Starting to set one up, writing the full description now.
Had to help Chris this afternoon, he was having more problems with configure.in
and Lapack on HPCC, but I believe I finally nailed this bug once and for all.
It now works for him.
Monday, Dec 3
Checking out mupad and PyMaple this morning...
Although I didn't immediately find a Python interface to mupad, I did find this guy:
http://swiginac.berlios.de/
Swigimac might be what we need, looking at documentation.
Swigimac also uses symbols, so I worry again as to what might happen with indexing vectors
Yeah, it's the same problem that sympy had.
I also found SAGE:
http://www.sagemath.org/
This may actually do the trick because it does not use Python string symbols.
It depends at a low level how things are computed.
Oh wait, yes it does. You have to define independent variables as objects of type var,
whose constructor accepts a string.
Well all I can do is try it and see... no examples that I can see in the tutorial.
News flash...takes forever to install. :p
Also, removed SystemFactory from MDL - put the build() functions for
Physical ad Forces with their repsective class definitions - there was reaally no need
for this factory - it actually wasn't a 'factory' really.
Now the factories which are left are Propagator and Force, these are both necessary.
Come to think of it, a helpful thing which I can add would be the ability to add a Python force to the factory.
And in reality, I think the force factory may even need a better design anyway. I don't use dictionaries like
I do in the PropagatorFactory, which leads to long conditional clauses. Will have to think about this some more.
Yes, I am going to attempt in the next couple of days a redesign of the ForceFactory. That will be necessary
before a publication. It's good that I'm writing this; I'm finding a lot more useful features for MDL.
Just did something pretty sweet- I got rid of prop.timestep() for MDL propagators. Now
you can just call self.dt. I create a binding for a member variable dt when the object is instantiated by the factory.
Sunday, Dec 2
Fixed compilation troubles with ProtoMol this morning:
1. Hessian w/ g++
2. Outdated Atom and AtomType classes
3. scpenergy-app
Also committed the standalone application to tools/.
One thing I did not do last week is commit the MDL with the updated Python force prototyping, so
I'll do that next.
Committed... The force prototype is in HDForce.py; which currently is in the MDL root since it's the only
prototyped force; I did not see the point of creating a separate directory for now although we'll need to eventually.
Simulation file is simulations/Impulse_Butane_HD.py
Also committed some of the libraries I used to compute energy gradients from Hinsen, they are open source.
Changed the README to specify the website I downloaded them from; and explained there were some minor
modifications.
Now have to work on my slides for class tomorrow...
Friday, Nov 30
Almost got the scpenergy-app done, there's just a few issues I'm having in setting force parameters,
such as cutoffs, etc. I'm trying to get around having to invoke the factories, but I may have no choice.
I do have an application running, it just produces zero for energy and force since the cutoff is currently zero.
Got the application working.
Working on the standalone per Dr. Izaguirre's request.
Close to achieving it, but not quite - few more linking errors but I have to go right now, will finish this weekend.
Wednesday, Nov 28
Looking at the issue from last night. Something's getting deallocated that shouldn't be.
Fixed, another dumb reference issue. Okay, it works now.
Beautiful. MDL propagators now work both as outer and inner propagators.
Now I'm going to return to Pmv because I know that interface is totally broken now, and I really
would like that module workingagain. Plus now that MDL is a bit easier to
use, hopefully my MDLCommands.py module will work better.
Newest version of Pmv works with Python2.4 - eek.
Python2.5 breaks the installation script.....
This wuold imply that if somebody is using MDL with Pmv, they['ll need Python2.4. And I wonder if we have back-compatilbitiy.
MDL breaks with Python2.4; and actually yes even Numpy 1.0 does not install properly with Python 2.4....
Okay they have different tarballs for tdifferent versions of Python. So I have to install the one for Python 2.4.
I'm going to see if I can run and install Numpy with Python2.4 and somehow then run MDL with Python 2.4.
Or would it be better to try to reconfigure Pmv to work with Python2.5? Now that I think about it, I remember Dr. Izaguirre could not install Python2.4
on his laptop without uninstalling 2.5.
If the tweaks are small, perhaps I could just make a supplementary tarball of
Pmv which works with Python 2.5.
Yes it was just one small tewak in an 'install' file. I'm just going to do that.
Okay, so it appears Pmv cannot run over an SSH connection. I tried that and my whole operating system crashed.
I'll try on LCLS... nope, Tcl and Tk headers are not installed. Oddly, the libs are.
Reinstalling the way it was, and I guess I'll just have to go work directly on a first-floor machine.
Tuesday, Nov 27
Merging the SCPISM Hessians to Protomol this morning...done.
Taking a little time this morning to go through some online documentation on PyFacebook...
Few class-related things for tomorrow, then MDL the rest of the day.
After fixing several more bugs, my HMC with RMT terminates with:
Step : 8, Time : 4.000 [fs], TE : 4.1359 [kcal/mol], T : 462.5085 [K], V : 121.71 [AA^3]
**** Move accepted
**** Move accepted
Step : 10, Time : 5.000 [fs], TE : 0.0000 [kcal/mol], T : 0.0000 [K], V : 121.71 [AA^3]
**** Move accepted
**** Move accepted
Step : 12, Time : 6.000 [fs], TE : 0.0000 [kcal/mol], T : 0.0000 [K], V : 121.71 [AA^3]
*** glibc detected *** corrupted double-linked list: 0x0a241508 ***
Abort
Monday, Nov 26
Busy in office hours (they have an assignment due at midnight), but I did fix the bug using MTS with methods.
Oops wait, nope - my energy is always zero - bad.:)
Fixed it - wasn't a bug, just was passing things in the wrong order to nextInteg() in the impulse() routine. Looking guten.
Committed the changes.
The other problem is this business with outer MDL propagators.
Sunday, Nov 25
One thing I'm realizeing this morning... the Propagator Factory
is not a singleton. It really should be.
Done! That was easy.
I think I also want to get rid of this SystemFactory class. It just contains two methods,
one to build the Physical and the other to build the Forrces, but they really can just be added to the appropriate classes I think.
I'm also fixing up the PropagatorFactory a little bit...I don't like the fact that
we use four
parallel arrays - I believe a Python dictionary should be used for the resgistry.
Also changing parameters to Python dinctionaries.
Oops!! We can't use dictionaries here. Because order matters when we pass parameters.
And Python dictionaries are like STL maps - keys are stored in sorted order.
NOW - we could compromise. Keep the registry as a dictionary, but use tuples for propagator
parameters.
Got it - okay that was a good comrpomise.
Man the factory is nicer now.
I'm thinking at this point pymat is not so nice anymore. You need to have Matlab installed to run that,
and Matlab is not open source. It is likely that
most MDL users will also have Matlab, but if MDL is going to be open source
I'd of course prefer any dependencies to be as well.
Switched to writing for a while this afternoon... I'm realizing it's going to be
delicate how I explain this PropagatorFactory. This easily could get confusing if I'm not careful.
I'm going to try to use some visual aids to help myself out on this one.
While writing this section, appears there are more bugs... for instance, I cannot do MTS on methods
anymore... because I'm now having create() return an application of the method,
the code no longer knows to call nextInteg() at the appropriate time. This needs to be fixed.
I'm not sure I can do this recursively, because the application of the next propagator is delayed until
nextInteg() is called frm within the method.
Also seems to no longer like outer MDL propagators, although I can use ProtoMol MTS with MDL STS (i.e. Impulse with RMT).
Mixing propagator objects and methods is not practical because typechecking is going to fail for SWIG-wrapped propagators,
for instance. The nextIntegrator data member, which is referenced all the time for MTS schemes, will expect an Integrator pointer
and not a method handle.
Meanwhile, I'm trying to set a minimal set of goals that I would like to get done by Christmas:
1. A good solid writeup on MDL that we can turn into a submission to JCP
2. Get the VMD and Pmv interfaces working again
3. Incorporate Matplotlib easier and more conveniently
4. (With a little help from above) the Python forces.
Next semester will be producing a Windows version and perfecting FTSM. I'm thinking we may
not even need these for a submission; since we already have complex schemes like NosePoicnare and RMT implemented in the language, also
Dr. Izaguirre did parallel-in-time tempering and Takahashi.
They would be nice things to have though...
Friday, Nov 23
Revisiting the mysteriosu force declinde for my Python harmonic dihedral force...
I have checked and the energy is the same, but the forces on each atom
are ever so slightly off, after a few digits. However, this difference is
clearly propagating.
It has to do with the method Hinsen used in calucating the gradiets
unfortunately it's hard to discern from the code itself and there
is not much documentation available.
I have delved deeper, the differeence in energies is definitely resulting from
the change in kinetic energy resulting from the velocities calculated by a slightly
approximated value for the forces. I'm afraid this library is no good
either.
What I really need is a cross between Hinsen's code and sympy: something like Matlab,
which can compute derivatives using symbols like sympy, but these symbols can be vectors like
Hinsen's code. Thus I could get a function back of say a vector x, then plug in my x as a vector.
Now there is this pymat, which is a Python interface to Matlab. Could this come in handy?
This also allows plotting in a Matlab window. We kind of have that with Matplotlib...
One thing I do worry about, PyMat is old, about 1999. The newest version I found
works with Python 2.2(!!) and Matlab 5. This may not work with Python 2.5.
Although ideally Python would be back-compatible.
It does have a sourceforge site, but sourceforge appears down for the holiday. :(
I went to the 'Cached' link in google, appears they have a 2003 version there, which works with Matlab 6.5 and Python 2.2.3
Trying the one I could get to, off a Professor's website...
No, these only work for Windoze. The sourceforge site has tarballs which can work with Linuxl
but I can't get to them at the moment.
Found another one on softpedia.org - followed the instructions for compiling on Linux
Yeah with Matlab 7.5, a LOT of these libraries are undefined - several linking errors and a lot
of symbols not found.
I was able to get them from different locations, but still cannot resolve the following:
/afs/nd.edu/i386_linux24/opt/und/matlab/7.4/bin/glnx86/libut.so: undefined reference to `__cxa_get_exception_ptr@CXXABI_1.3.1'
Guess I'll write a little more, maybe I'll get lucky and sourceforge will be back up soon.
Wednesday, Nov. 21
Have fixed several bugs in my force prototype but still the energy is mysteriously declining when
I used my Python harmonic dihedral force versus the SWIG-wrapped version.
Somehow, there is an empyt hook I think.
Anyway, taking a little time to write this afternoon.
Tuesday, Nov 20
I've found something interesting:
http://ccl.net/cca/software/SOURCES/PYTHON/HINSEN/README
I don't know what this software is, but I wonder if we can interface to it.
Seems like its purpose is for computational chemistry anyway. Looks like
it's called Hinsen or that's the uauthor...
Nope, this gives me the same problem. It doesn't like that I index a vector
that I'm trying to differentiate with respect to. :(
I like these libraries better than sympy however.
*May* be able to pull this off with four independent variables.
Okay I've almost got this working... the biggest issue I'm having right now
is the type conversion between Hinsen's vectors and numpy arrays (analosogous to the problem
before with Vector3DBlocks and numpy arrays.
Almost got something set up but my energy is not correct. I like these Hinsen
libraries however.
Monday, Nov 19
Trying to work out the derivative this morning. One issue is that
the dihedral function is actually not a continuous function of x, but piecewise. Sympy's derivatives substitue independent variables by symbols, so I can't use x
as a vector in this function.
Yep, this differentiation does not work if the dependent variable is a vector and we access certain
portions of the vector. In the same way, even with respect to coordinates of individual atoms,
from a computational perspective these too are vectors. And when I compute dihedrals, I have to do cross products.
So theoretically this is a function of twelve dependent variables (4 atoms, 3 coordinates). Ugh.
Exploring for better options before giving up...
I have tried passing symbols (i.e. x1, y1, z1) as arrays, but there are two problems:
1. I cannot pass them as a Python array because operator - (for instance) is undefined.
2. I cannot pass them as a numpy array, because these arrays must hold numbers, and the moment I replace them
with a symbol I get a type mismatch.
So it's definitely the case with sympy derivatives, that *any* dependent variables, whether they are contained within
a vector or not, must be specified as symbols.
Now, what I may be able to do, is add this behind the scenes. Meaning, from the user's perspective,
they provide the energy function as a function of the positions x. They then call my routine for differentiation.
My routine, behind the scenes, performs some kind of translation, obtaining each dependent variable and
performing the differentiation, translating back into a force vector.
This is without a doubt going to be challenging and I have to think about methods for designing this.
Saturday, Nov 17
Fixed the following:
1. Added += and -= operator overloads to the Vector3DBlock class which accept both
Vector3DBlock and Vector3D.
2., Changed restoreForces() and saveForces() from shallow to deep copies of data
3. Added an MDL interface to XYZTrajectoryWriter, passing a 1 to the tuple
in the io dictionary will print the old forces instead of the current forces.
4. Fixed another shallow copy in the OutputCache
Thursday, Nov 15
Working on this OutputXYZTrajectoryForce modification this morning...
Finished.
Now going back to the forces... I actually got something working with
this approx_fprime but unfortunately, the numbers seem inaccurate. Total energy goes
down - not extremely fast, but fast enough to know something's wrong.
Oh, I just needed to use a higher order approximation.
Okay, with a 24th order approximation I was able to obtain forces
that were relaitvely close (within a few decimal places), but still I'm not sure that' going
to cut it honestly.
It would really be ideal if I could find a Python library devoated to diffeerentiation and not Scipy which is humongous and kind
of a jack of all trades but a master of none.
SymPy might be an option...
Okay, here I've got something:
>>> def y(x):
... return x*x
...
>>> def z(x):
... return x*y(x)
...
>>> x = Symbol('x')
>>> diff(z(x), x)
3*x**2
So if I can just find a way to apply the final expression with the symbol x replaced by some value...
Okay, yes - I can convert the expression returned by diff() to a string and then call eval:
>>> h = diff(z(x), x)
>>> eval(str(h), {'x':4})
48
>>> eval(str(h), {'x':3})
27
So this is good -- gonna try this now, will be better than just
making an approximation.
Wednesday, Nov 14
Okay, so numpy has a function polyder() which can compute the derivative of a polynomial,
given an array of coefficients. Obviously this doesn't help too much, we'll have more complex energy functions
and have to execute the chain rule, etc.
I think scipy may have a better one however. I'm going to look into this.
Tried installing, but it appears scipy is not compatible with the intel fortran compiler:
Fortran f77 compiler: ifort -72 -w90 -w95 -KPIC -cm -O3 -unroll -tpp5 -xM -arch SSE2
Fortran f90 compiler: ifort -FR -KPIC -cm -O3 -unroll -tpp5 -xM -arch SSE2
Fortran fix compiler: ifort -FI -KPIC -cm -O3 -unroll -tpp5 -xM -arch SSE2
creating build/temp.linux-i686-2.5
creating build/temp.linux-i686-2.5/scipy
creating build/temp.linux-i686-2.5/scipy/fftpack
creating build/temp.linux-i686-2.5/scipy/fftpack/dfftpack
compile options: '-c'
ifort:f77: scipy/fftpack/dfftpack/dcosqi.f
ifort: Command line warning: extension 'M' not supported ignored in option '-x'
ifort: Command line error: Unrecognized keyword 'SSE2' for option '-arch'
ifort: Command line warning: extension 'M' not supported ignored in option '-x'
ifort: Command line error: Unrecognized keyword 'SSE2' for option '-arch'
error: Command "ifort -72 -w90 -w95 -KPIC -cm -O3 -unroll -tpp5 -xM -arch SSE2 -c -c scipy/fftpack/dfftpack/dcosqi.f -o build/temp.linux-i686-2.5/scipy/fftpack/dfftpack/dcosqi.o" failed with exit status 1
At least not the version we have installed.
I may try a lslightly older version just to see.
Nope. The old one doesn't work with Python 2.5. :-)
Okay, I installed Intel Fortran version 10 and we worked. It's only a 30-day evaluation version, but fortunately now
I should have Scipy permanently.
Okay yes the one in scipy is a bit nicer; this can compute derivatives of functions of x:
>>> from scipy.misc.common import *
>>> def y(x,p):
... return p[0] + p[1]*x + p[2]*x**2
>>> derivative(y, 2, args=([[0,0,1]]))
4.0
But unfortunately it produces unexpected behavior if I have other functions of x in the equation. This should
be 12 for instance but for somre reason it is returning 13:
>>> def z(x, p):
... return x*y(x,p)
...
>>> derivative(z, 2, args=([[0,0,1]]))
13.0
There may be a way around this I still have to look more.
I do see some gradient functions available.
There is a function approx_fprime(). Depending on how accurate this is,
it could be a good function to use.
Looks like there's also an approx_hessian() function...
Monday, Nov 12
Yesterday was bogged down with stuff for class, today I'm going to produce a prototype of the
Harmonic dihedral system force using what I currently have with Python. Again, this will not currently
be pretty, but I would like to figure out what kinds of areas I can simplify.
Ouch, also this is not too powerful at the moment - there is no way to take into account
boundary conditions (without hardcoding them anyway). I don't really want to call minimalDifference - ti would be nice to just be able to do
for example pt2 - pt1.
Alright I have written up an HDForce, and it calls the valuate method and updates the forces,
but somehow my results vary from using the wrapped HarmDihedralForce.
They are very close though... to the point where I think it's actually an issue of precision..
I'm leaning this route because at the first step things are identical to eight digits, but after the first about the eighth
digit starts changing. Then it 'propagates' its way up.
Saturday, Nov 10
Okay, I figured out the problem with the PyDict_GetString() routine. Member functions are not
actually stored in the dictionary of an object, only member variables. Which explains why "evaluate" came back NULL.
So I need to figure out how to access member functions of an object.
Member functions are however in the dir().
YES!!!!!! Alright!! I was just able to call a member function on an object using PyObject_GetAttrString().
I added a new Python force object which simply incremented a counter in its evaluate() routine, added it
to a back and force list in ProtoMol's ForceGroup class, and the evaluate function was called repeatedly and held
the values of the member variable counter, because I rpinted them in each call to the Python evaluate method.
And after the propagate() call the final value of counter was correct.
So this means... we probably can add forces now!!! Alright exciting!!
Now, the downside is -- currently, we'd pretty much have to set forces and energies the same way as with ProtoMol.
We'd have to accumulate the force at every point in the atomic force vector, then compute the energies sepeartely.
BUT, we now have Numpy. So I wonder, if there is a function which can compute the gradient of a function
If so, we could get away with only specifying the Hamiltonian. When we take the gradient of the Hamiltonian with respect to the
atomic postiion x, we should be able to just plug and chug to get force values, and we should be able to do it all behind the scenes.
Yes - absolutely, even in the case of a harmoinc dihedral force in ProtoMol, that could be accomplished with one statement for the energy,
but the force on each atom takes a lot of work (about 95% of the function).
The tough part is, if there are other functions involved which depend on the position. For example, in something
as simple as harmonic dihedral, even the dihedral angle phi is a function of x. So getting the gradient of the potential
becomes a matter of apllying chain rule: dV/dx = dV/dphi * dphi/dx. If there are more functions of x involved, we get even
more complicated. So hopefully, numpy would bea ble to deal with this.
Anyway, off to the game...
Friday, Nov 9
Okay, I'm at a bit of a crossroads with these Python forces...
Here's the issue. I can in a .i file for SWIG, extend a C++ class and define routines to accept
Python objects and call their functions. However, these routines are only invokable
through Python (they are not recognized when the class implementation .cpp file is compiled).
Therefore in a nutshell I cannot call a routine that I specified in a .i file from a .cpp file because
the compiler does not see the function (it's not in a .h file that has been included).
Now in the .cpp file, I can define routine(s) there which operate on Python objects, but it appears
in that situation a Python interpreter must be initialized. So we have two running. :-(
I want to have both. I want to be able to, say, from a back end calculateForces() routine, invoke a routine
in the .i file (thus not haveing to start a new interpreter) which operates on Python objects. And at the moment it is not clear to me how.
Although I get the same error (seg fault, Null function) in the .i file.
Looking at Dr. Skeel's notes again..
ALRIGHT!!! I got a simple Python function call to work from C++.
And I was using the back end to call it. So definitely, I don't need to start a new interpreter. Phew. What a relief.
So now, this is an unbound, non member function with no arguments.
I also can pass primitive type arguments.
So the next step now, is to be able to call a function of an object.
The thing is, I will not be able to, from the C++ backe end, pass non-primitive data
to the Python routine. When setting up the tuple for parameters, it requires primitive data types.
BUT, I can just have a force object from Python contain references to
the position, velocity and force vectors. That way any changes made
by Python will be seen in the back end, and vice-cersa. So my meber function
does not have to take parameters.
In fact a new force object can just contain Physical and Forces references. Bene.
So, returning..how do I call a function of an object.
When I try this:
PyDict_GetItemString((*currentForce), "evaluate"), where currentForce is a PyObject*
representing an object of a Python force class, this returns a NULL function.
Thursday, Nov 8
Took care of a few small things this morning, now working at the dynamic allocation of SCPISM
parameters....
The way I see it, we have two options: I could (1) make a new keyword for ProtoMol 'doscpism',
which if set to true builds SCPISM parameters, or (2) I could build things in the construction of SCPISM
forces. I really like the first option a bit better, because I could be using a CoulombSCPISMForce
or CoulombBornRadiiForce, and each one needs the parameters built. But I don't want to do it twice, if I include both forces.
Completed late this afternoon...
Wednesday, Nov 7
Santanu got an error with the new ProtoMol (with SCPISM), but I was unable to reproduce it
after a fresh checkout.
Oh but I was with g++. Fixed it.
Wrking on the force hooks, currently getting a seg fault... which is because the evaluate() routine
that I am defining is for some reason coming back NULL.
I had to take some time to grade this afternoon, so I'll continue this tonight or tomorrow.
Tuesday, Nov 6
Made a copy() method for physical which dose element by element and committed it...
Now my focus is going to be these force hooks. This is something missing right now.
So, the current issue right now is that I need to be able to add new force objects to an MDL
ForceField (which inherits from ForceGroup), define their evaluate() methods and make them callable
from calculateForces() invocations.
Clearly, I'll hvae to at least override calculateForces(). I'm looking at the SWIG manual; again to see if I can
come up with some creative way to call Python functions from C++ in the .i file.
Gasp! I believe I just got the ISM code working in ProtoMol.
Validating...
Yup! Alright!!!!
Committed and let the group know...
Monday, Nov 5
Comitted the new IO interface...
Another thing I want to fix once and for all is this memory leak when propagators change..
Gasp! I finally fixed it - wow that was tough.
Off to dinner... little grading to do tonight.
Setting off another run with kappa=15.
Sunday, Nov 4
This morning, cleaning up the interface for the forces, using a similar method with the dictionary.
Great - got them working! They look excellent!
I also have the dirty bit working, such that forces are only rebuilt as necessary.
And I removed addAllForces(). That gets rid of an extra level of overhead.
Committed the code, my next goal will be to clean up the IO a bit.
Saturday, Nov 3
Fixing up MDL, I think I have found a way around passing phys, forces, etc.
every call to propagate. I could just pass them to the Propagator() constructor. Then all people will
have to pass to propagate() is the method, number of steps, timestep, and forcefield. I'll need to keep
forcefield there however, because if there is an MTS scheme.
Here's what I have working this morning:
prop = Propagator(phys, forces, io)
prop.propagate(scheme="BBK", steps=200, dt=0.5, forcefield=ff)
Which seems much nicer. Nwo I have to think about MTS cases, where theoretically there can be arbitrarily many:
What I've previously done, is allowed the user to specify additional arguments, i.e.:
prop.propagate("BSplineMOLLY", 200, 5, ff, "Leapfrog", 0.1, ff2) # cyclelength=5, dt=0.1
What I can do now, is accept lists for the scheme, cyclelength, and forcefield, so it will hopefully look like this:
prop.propagate(scheme=["BSplineMOLLY", "Leapfrog"], cyclelength=5, dt=0.1, forcefield=[ff, ff2])
The rules will be:
- Each element that is passable as a list should also be passable as a single value as appropriate (i.e. for STS,
I should still be able to do forcefield=ff for instance and not [ff].
- The way to know the number of levels will be by the cyclelength. The number of cyclelengths provided equals
the levels of propagation - 1.
- If there are fewer schemes than cyclelength, we should assume the innermost propagator to be Leapfrog and
outer propagators to be Impulse.
I like this a great deal. I'm going to try to get this working...
The only thing I also have to think about: extra arguments. Where should they be supplied?
I've got it! I could use Python dictionaries, like so for normal modes:
gamma = prop.propagate(scheme=["NormModeInt", "NormModeMin"],
steps=1250,
cyclelength=1,
dt=4.0,
forcefields=[ff, ff2],
params={'NormModeInt':{'fixmodes': 40, 'gamma': 91, 'fdof': 0},
'NormModeMin':{'avModeMass': 30}})
This looks nice.
Got it working!
Okay, now I'm getting more ideas with this dictionary.
Why not, use it to set data for LennardJones and Coulomb? This would save us SO many lines
calling accessor routines of ForceField.
Even with the output we could manage it. Something like:
io.plots = {'total energy': 4,
'temperature': 2}
Anyway, I'll finish with the propagators first. I like this so much better than calling all these methods though.
Looking at the PIT, it seems that I'm going to have to be delicate in terms of when I rebuid these Physical objects.
I think the rule should be: associate a dirty bit with every Physical object, initially true.
In a call to propagate(), if the dirty bit is set, rebuild the Physical object and clear the dirty bit.
Otherwise don't rebuild. tiem will be the one attribute exception that does not set the dirty bit.
Anyway, game time - will come back to this later or tomorrow.
Friday, Nov 2
With 20 didn't buy us much:
Thinking, I could change it to 10, or I wonder if I have to lessen it even?
Thursday, Nov 1
Okay turning up kappa definitely seemed to make the string worse....
It's actually unclear though if it really made it worse.... so I increased kappa by a factor
of 10 - what does that mean? So it means the constraint is 10 times as strong.
But it also, on the flipside, should provide a more significant update of free space.
Maybe this constraint is just too strong though.
Let me go to 20.
While that runs, having another look at the MDL code..
One thing I'm not a big fan of is that when I create a force object, I have to add it both to
a Python array and to the back end structure.
Ideally I could just add it to the back end. I wonder if ForceGroup
has someway to access individual Force objects so I can change their parameters dynamically.
One thing I can do is get rid of this DihedralIIISystemForce, since that was used for old FTSm
AH! Okay, I think I found a way to cut the size of MDL significantly!
I think we may only need the shared objects and not the Python modules that get created when I use
SWIG.
That is going to remove a lot of confusion I think.
Well I can't do it in all cases... if I define a template as a class, that is undefined in the shared object binary.
Anyway, this afternoon I have to change gears, because I need to present my C++ class with a studyguide for tomorrow.
And I have to do this away from the office, because the windows box is now down.
Wednesday, Oct 31
Office hours completely flooded this afternoon.. I did manage to throw in some more pydoc into
certain classes and make the factory interface a little more clear in MDL .. not committing anything til
it's fully functional though.
Tuesday, Oct 30
No good on the string....
I still think everything is correct, but I may need to bump up gamma periodically like Luca did.
2000 for gamma is clearly too small.
Turning gamma up... 20000 is still too small.
This seems way too high, maybe something else is wrong sigh...
Working on some of the other suggestions on the MDL Wiki page..
Committed a new Propagator.py with a propagate() function that does not return positions and velocities.
Added a query() method to the PropagatorFactory.
Committed the new directory structure for MDL.
Monday, Oct 29
My fears were confirmed with FTSM, we definitely have a memory leak someplace.
I had suspected this when I saw how slow it was running, and especially
since it was getting slower with time..
Traceback (most recent call last):
File "simulations/FTSM_ALA.py", line 109, in <module>
M = FTSM.M(physarray[workpt][0], PHI_DIHEDRAL, PSI_DIHEDRAL)#, physarray[workpt][0].phi(PHI_DIHEDRAL)-S[workpt][0], physarray[workpt][0].phi(PSI_DIHEDRAL)-S[workpt][1], kappa)
File "/afs/nd.edu/user25/tcickovs/Research/mdl/FTSM.py", line 136, in M
rij = phys.positions[atomJ*3:atomJ*3+3] - phys.positions[atomI*3:atomI*3+3]
MemoryError
We only got about 46 ps with unsolvated alanine, which I don't expect to be anywhere near long enough:
As expected, the string has not moved much although it does appear to be going in the right
direction.
Trying a new method where I do not redefine new HarmDihedral objects but call a member function
toset the reference dihedral values.
While that runs, I'm looking more at this ISM merging.
I feel convinced at this point, that we are not going to be able to get around
making the data member topology of OneAtomPair a non-const.
We pass it to every nonbondedForceFunction including CoulombBornRadii
force, which will change this object. The only other option we have is to move the bornRadius, sasaFrac, etc. out of the topology.
Which may not be an insanely bad idea actually.
We've definitely still got the leak in FTSM. It's running horribly slow even after 16 ps.
Looking at top, we're definitely increasing by about a meg of RAM every couple of seconds. No good.
Probably what's happening is that eventually swapping is necessary when all the RAM is taken up,
and the swap rate is simply too large.
Okay without a doubt the propagation of x and y is causing the leak.
Without it, we hover right around 140 MB.
Put in a quick fix for it, the leak was in creating new Langevin propagators for each replica at each
step, I have now hardcoded to only do it the first time propagate() is called.
This is not a permanent solution, but I'm more concerned with getting results so I'm going to note
the hack and now look for a better way.
I think the BEST way would be to actually keep a cache of the current propagation scheme, timestep, forcefield,
variables, etc. Then only create a new propagator object IF some of these change.
I also have to figure out though, why we have the leak so when I do the above we don't have issues.
Friday, Oct 26
I created the three types of forces this morning.
There is still an unforseen hurdle though, OneAtomPair actually has
a const data member for the topology. It then passes this data member
to the operator() function of whatever force is being evaluated.
This is a bigtime ouch. Although I can define two functions one of which
takes const and the other which tkes non-const, doing this to a data member is a little more
challenging....
Now, I could just make it non-const, but then *every* wrapped force would
have to have an operator() which took a non-const, even if it invoked the const. That's a mess.
Anyway, in order to free up Monday, I have to take care of some class stuff right now.
This week I have to do the following:
1. Make a new homework
2. Make the slides.
3. Make a study guide for exam 2.
4. Make exam 2.
Busy busy...
Thursday, Oct 25
Looking at my new FTSM results this morning, I get a singular matrix crash
after 10 ps, again when I started bumping up gamma.
I'm starting to not believe that we can assume anything
about these kappa and gamma values.
Luca and Eric ran on a solvated alanine dipeptide system.
This afternoon I am working a bit with ismprotomol, and am going
to attempt to get rid of that const problem (my book came yesterday).
Okay, it is clear to me how this pattern works.
Taking another look at ismprotomol...
Okay I see how we can incorporate this change. We will however have
another small issue that I failed to consider... we defined a NonbondedCutoffBornForce
to compute born radii. So I assumed by doing this that we would be using a nondirect
algorithm for pairwise evaluation. Now that assumption is down the tubes.
I of course do not want to duplicate code and have a NonbondedSimpleFullBornForce
and a NonbondedFullBornForce.
Merged in CoulombSCPISMForce - that was the easy part.
Now with CoulombBornRadiiForce, that is ALWAYS right now paired to a NonbondedCutoffBornForce.
Theoretically what we want, is that you can simply slide this in to any template and have it work.
But now, the templeate has to potentially run a routine which updates the forces and born radii.
I have to do this *very* carefully, and must make sure that any portion which invokes this
code snippet only operates on a nonbondedForceFunction of type CoulombBornRadiiForce.
Otherwise, everything will crash (dR won't be defined, etc.)
I think we have no other choice here but to create NonbondedCutoffBorn, NonbondedSimpleFullBorn
and NonbondedFullBorn. Otherwise, we don't have any idea what TForce gets initialized to,
and cannot invoke any kind of functionality to compute Born radii.
I'll look at this again later.
On another note, I increased gamma for the whole situation and now have about 32 ps of an
uncollapsing string, but the simulation is starting to go very slow.
This is interestingly even the case on stones, 64-bit and I don't think there's any kind of process management
and decreasing nice values depending on runtime.
I'll give LCLS a shot too..
Ah, the good old libstdc++.so.6 problem on this box... fun times.
I really ought to be running in parallel on ISS anyway.
Although I wanted to get serial working, it seems that it is just too slow to get results in any
reasonable amount of time.
So that'll have to be one of my next goals.
Wednesday, Oct 24
41 ps this morning; still no collapse but I expected it to run for longer than that. :p
Yeah clearly not long enough, still the string hasn't moved much.
Something is wrong, there appears to be a memory leak. I stopped it after 42 ps,
we had 34 ps when I left yesterday, so clearly we are getting slower with time.
That should not be happening.
Okay this was my own stupidity .. I had turned off all other forces because I was
testing PM's Harmonic Dihedral force.
BUT - I found another bug! When I was computing M, my dihedral indxing was off by 1!
Fixed this, and we're looking better in the early going, see if it holds up.
Okay, I have confirmed now that the derivatives I am computing are the same as
the harmonic dihedral force. Fixed a bug (my atom indexing was off by one).
Now it's running and it looks mucho mejor (much better)!! I'll let it go.
Nope, still collapsing. :(
Gonna try and bump up gamma.
Tuesday, Oct 23
Santanu and I are working this morning on FTSM stuff...
We have a two-way method of testing ... I'm going to try just simply using
the Harm Dihedral force in my simulations and making sure that results correspond.
Man I've got a mess.
Harm Dihedral and my constraint each give me different results.
However, both give me a string which increases positive definitely in terms of phi and psi.
So each are not alike, and both give me wrong results.
Looking at my M computation...
Ah, okay! I fixed a bug in the constraint - I'm using harm dihedral now, and I was
accidentally setting the same target phi,psi on each replica.
Letting it run, and after 4 ps still no collapse. I will plot after 200 ps like Dr. VDE did.
Friday, Oct 19
I have an incredibly simplified Makefile.swig and fullCopy.sh, but Duong is having
some problems so I'm helping him out..
Gasp, finally finished them. They'll be easier than b4.
Thursday, Oct 18
Makefile.swig is much simpler per the early afternoon, going to give it one more compltet pass
to make sure I did not break anything.
Also, at some point I'm going to have to revalidate the VMD and Pmv interfaces....
Tried to install VMD from source and it refuses, it tries to include some header files which
don't exist. :( The binaries I have do not include Python support.
Wednesday, Oct 17
First thing this morning, is to set up this Wiki page.
Done here.
Now I want to brainstorm a little bit about how to better automate the Makefile.swig/Makefile.am
process.
There are a few things I see here:
1. The SWIGPYTHON #define is sometimes used to determine code which is necessary
for SWIG-wrapping but not when running straight-up ProtoMol. So actually problems could come in
if a user has just compiled linprotomol and then simply runs 'make -f Makefile.swig'. If the .o files
already exist, they will not be remade with the changed #defines. Sadly there is no way around the #defines,
because there is certain code that is simply not wrappable.
2. It seems that I am using the same command to compile all object files. I remember with BioLogo
there was definitely some way to automate this. Also I remember Maciek doing something with CC3D, automating
swig wrapping. Going to have a look at that as well.
Although, he used cmake. However there was an old version with automake+friends.
I've cleaned up the object files in there...will shrink the file by quite a bit.
Now looking at automake...
It seems that the dependencies are the same for all integrators, all forces, etc. (the obvious exception
is the .cpp file containing the integrator we're working with.
So I definitely should be able to macroize these and make the file simpler.
One I have my Makefile.swig as simple as possible, then I can work on perhaps generating it through
a Perl script or something.
Tuesday, Oct 16
Tried to set up an MDL Wiki page this morning, but it was not clear from
the website how to do so. E-mailed their webmaster.
I ordered Meyers' book last night, which should be arriving this week.
Now that we have eclipse working, I am more and more liking the idea
of merging the linprotomol and mdl repositories (well, more like - make linprotomol
a part of the MDL repository). That way both can be developed.
A disadvantage is that one cannot compile the C++ code from the Eclipse IDE,
though one can develop there. It will still be necessary to run either Makefile.swig
(for wrapping) or automake+friends, etc. for a standalone executable.
I've discovered one problem with the BBK run ... initial velocities are different somehow.
Even though we have the same temperature and seed...
Okay, I've figured it out. Even though I initialize two separate velocity vectors for each
Physical object, the back end only uses one random number generator.
Alright, I fixed that problem and now have idential velocities to start, but I still get different final positions.
Okay, well it's the same deal. The two of them use the same random number generator, so after you run one fora while it's off again.
If I could just somehow restart that generator. I think I may to thrive off the checkpointing.
Got it, I just had to reseed the generator when the integrators were intiializaed.
Yay!
Okay, now on to bbk().
Got one iteration but for some reason it seg faults later.Could again be an issue of copying.
Yep - that was it. It runs now.
Recommitting... had to debug randomForces() slightly and move it to the Forces module.
Also, fixed Propagator.propagate() to return a tuple of copies of the initial positions and velocities.
Monday, Oct 15
Getting an error when installing eclipse:
Warning: -Xms40m not understood. Ignoring.
Warning: -Xmx256m not understood. Ignoring.
Warning: -jar not understood. Ignoring.
Exception in thread "main" java.lang.NoClassDefFoundError: .afs.nd.edu.user25.tcickovs.Research.eclipse.plugins.org.eclipse.equinox.launcher_1.0.1.R33x_v20070828.jar
at gnu.gcj.runtime.FirstThread.run() (/usr/lib/libgcj.so.5.0.0)
at _Jv_ThreadRun(java.lang.Thread) (/usr/lib/libgcj.so.5.0.0)
at _Jv_RunMain(java.lang.Class, byte const, int, byte const, boolean) (/usr/lib/libgcj.so.5.0.0)
at __gcj_personality_v0 (/afs/nd.edu/user25/tcickovs/Research/eclipse/java.version=1.4.2)
at __libc_start_main (/lib/tls/libc-2.3.4.so)
at _Jv_RegisterClasses (/afs/nd.edu/user25/tcickovs/Research/eclipse/java.version=1.4.2)
It seems like this tarball already has the binary there.
I'm wondering if I can rebuild it from source. Maybe Ineed a more recent version of java too,
the one on these clusters is old.
Ah, okay - I need the Java runtime environment.
Got it. Eclipse works, now I've got to try to figure out how to get all three in place.
Appears there may be an online tutorial on this...
Alright! I've finally got eclipse running some Python scripts.
Now, I've run into the first hitch with MDL .. it appears that scripts are run
using absolute paths and that the cwd for the IDE changes even if you run it from the MDL
root. This is problematic, because when I register the propagators with the factory, I look in the
propagators/examples directory using a relative path.
Although this probably ought to be absolute anyway. Actually, I know what I can do.
I should probably change the MDLSetup scripts to set an environment variable MDLROOT.
Then I can from Python get that environment variable and append it accordingly.
Ok, doing this -- there is another challenge I see though.
It appears that with Eclipse, it is necessary to put in the PYTHONPATH yourself.
This is not good. A typical user of MDL will not know the folders to put in the path.
That was the whole point of MDLSetup.tcsh. Now if I could somehow get Eclipse's path
to reference the path in the environment variable, we're in business.
New problem now:
terminate called after throwing an instance of 'std::bad_alloc'
what(): St9bad_alloc
Which unfortunately is behind the scenes. :(Well I see the problem, my position and velocity vectors are empty.
Ah okay...similar issue. Input file paths are relative.
So I should have a check at the beginning of my functions for file input,
to check if the file as it stands exists; if not append the MDLROOT.
YES!! Okay, I got a simulation working with eclipse, pydev and MDL.
Yep my MDL simulations run. This is indeed pretty neat. I like the hyperlink
stuff too.
I wonder if I can run IPython.
No, ipython is giving me problems. I cannot manually add the executable,
and when I add it through a file manager it somehow prints this long message "Successful update! ......"
which is not the name of the executable or course. :p
This IDE is pretty nice, the only thing I don't like is this PYTHONPATH
business. Unless there is a way that I don't see yet to automate its settings.
Anyway, okay this is good ... now my next step is going to be to check out the book
that Dr. Izaguirre recommended last Friday, to try to resolve this const/non-const issue.
The book is called Effective C++, Third Edition, by Scott Meyers.
Looks like I'll have to get it on Amazon, University has an older edition.
Will purchase it later; I'm going to create an MDL Wiki page on SimTK now.
Dr. Izaguirre, here are the instructions for running an MDL simulation with eclipse:
1. Start Eclipse.
2. Go to Window->Open Perspective->Other... and choose Resources. I've found that's the easiest view to use.
3. Go to Window->Preferences->PyDev->Interpreter-Python
4. Add the executable for python 2.5. Also, set the System PYTHONPATH. YOu have to do
it folder-by-folder using the New Folder button (I was looking for a better way - see above - but
have not been able to find one yet). My suggestion is to just source MDLSetup.tcsh (or bash), run
'printenv PYTHONPATH' and then run through the folders and add them.
5. Go to File->New->Project->PyDev->PyDev project. Give it a name (I used MDL), you can use
the default contents, make sure it runs python2.5, and leave the box for the src folder checked.
6. Right click the 'src' folder in the Navigator window, and click 'New->File'
7. Click the Advanced button in the resulting dialogue box, check 'Link to file in the file system', and then
browse for an MDL simulation file. When you find it, double click. It should now appear in the Navigator
window under the src folder.
8. Go to the Navigator window, right click on the simulation file, and click Run As-> Python Run. Should work,
with output displayed in the console window.
Let me know if you have troubles.
Sunday, Oct 14
Installing eclipse, pydev and ipython this morning...
I'll have to read how to work these things...
Printed up the ipython manual..
Friday, Oct 12
Very close this morning, my only error is at the very end:
*** glibc detected *** double free or corruption (!prev): 0x08f9ce88 ***
Checking this out. Once everything is functional as far as a standalone
executable, I'll retry copying things over to MDL and make sure I did
not break anything there.
Okay, got it. Now testing with MDL...
Appears MDL is still in place, running all simulations just to make sure.
Yep, we're good. Recommitted linprotomol.
Now, looking at the PIT* files that Dr. Izaguirre committed,
which have the crash at the end.
Also, I need to make a note: when I get time, install
eclipse, pydev and ipython for a graphical IDE and debugger.
And create a simtk Wiki page for desired features, bug fixes, etc.
Looking more at the PIT, it is clear to me that this crash comes
the moment we start propagating using two more more NormalMode
propagators. If I kick out just after propagating with the first one,
the program runs fine.
Yep, sure enough - the same thing happes in my basic NM runs, if I
add a second propagate() statement with the same NromModeInt,
get the same double free.
But it's not a general propagator issue, because I definitely can run BBK back-to-back,
or Langevin.
Got it! It was an issue of pointer copy as expected, we were doing
an *ex0 = *myPositions in NormModeInt. Changed it to ex0->intoAssign(myPositions).
Recommitted.
Thursday, Oct 11
Quote of the day: "Vector3D is like a cockroach that cannot be swatted".
Ok, I almost have it... it's down to linking errors finally.
Also helping Duong get linprotomol compiling...
I have my linprotomol application compiled now, but when I run it with any simulations,
it fails with an empty:
Fatal Error:
but no message. :(
K, IntegratorFactory::make flunked, but interestingly no error message was written
(the error message was what was supposed to be output here, but it's an empty string).
I've checked and the integrator string is correct.
Wednesday, Oct 10
Gasp... okay, now that I've finished with coursework I can *finally* put some serious time
into getting this linprotomol executable.
Debugged a lot of the new integrators tonight.
Tuesday, Oct 9
Okay, I have 26 ps before terminating with:
terminate called after throwing an instance of 'std::bad_alloc'
what(): St9bad_alloc
However my string is absolutely horrific. It is slightly collapsed,
and is a concave arch above the (originally) diagonal string and should be convex, below
the diagonal string. Bugs bugs bugs...
I'm going to check first and see if my Langevin propagator is accurately sampling phase space
at the two ends.
Ah... okay I may have found part of the problem.
It seems that C5 axial samples very well, but if I start a 10 ps simulation from C7 axial, the molecule
is dragged (and eventually stays) in the C5 axial. This is problematic. C7 axial should be a metastable state,
even if C5 is more energetically favorable. I have no constraints in the simulation right now.
So this needs some debugging.
Couple emergency items that I have to work on after making my slides:
1. Check out the PIT_... file
2. Make linprotomol create an executable again (it is broken right now).
Monday, Oct 8
Committing the code that Dr. Izaguirre asked for...
OH! It finally occurred to me why things were becoming unstable
after 10 ps.
That is the point where I start increasing kappa. So maybe for this system,
it's not a good idea to increase it like Dr. V did.
Trying this again...
Okay, now it runs 12.6 ps and becomes unstable.
Friday, Oct 5
I have fixed the reparameterization to allow the last point in the string to move.
Unfortunately when I do this, the entire string collapses into the first metastable state. :-(
i.e., it moves out of the second one.
I can tell this is because free space is moving too fast; that there is not enough time for
each individual replica to get into its appropriate state. We have to drag things across
a huge area of energetically unfavorable conformations.
But then again, it's still getting pulled towards the first. :p
Actually I solved this (at least temporarily) by taking a very large gamma (so that free
space is updated very little while Cartesian space has a change to grow into the right places.
But, it seems we may have a memory issue of some kind after 10.1 ps, my string is somehow
reduced to a length of one point.
Trying with fewer replicas...
It's going. I'm going to have to wait a little bit for it, so I'll switch to another machine and
keep doing the epydoc comments.
Okay, I see why the string is going to one point - basically the system becomes unstable
at around 10 ps.
Translation is probably that one of the replicas gets unstable and throws the rest off
in reparameterization.
But probably, if I could do a reasonable sweep on kappa, I could find values that
yield more stable dynamics.
Shrinking kappa for now...
Thursday, Oct 4
Opening the morning by making slides for tomorrow.
Done, now switching back to the FTSM, trying to run on the simple
alanine dipeptide case...
WOW. OK, THERE IS A SERIOUS BUG SOMEPLACE.
So I'm printing out the values of the string. For some reason, when
I switched from CoulombDiElec with a cutoff to straight-up Columb with no
cutoff, I get identical values in the string after six iterations!¡ So something
is seriously wrong in the control flow. This could be the result of something not being
initialized properly. I'm not sure at this point, but I plan on figuring this out today.
The force vector has correct values.
Fixed it! Was a bug in Langevin Impulse, because of the way I wrapped things the half
kicks were not being called.
Okay, now it's better, but still the string is collapsing.
I've checked the values for M and the collective variable difference.
The bigger problem is clearly the collective variable difference (M values do not increase
with time but the collective difference does).
Fixed that bug, was a +/- issue.
Now, the string looks a lot better in the arly going but crashes after 1.7 ps
after a failed reparametrization.
I am currently running with 40 replicas and I think that's too many (the points
are crossing). Going to 20...
MUCH better. The string actually is looking good (although I ahve to run for
longer). There is only one thing at this point that I'm unhappy with, and that's
the reparameterization routine -- it's fixing the right endpoint, which is wrong.
Goin to make an attempt at fixing that.
Wednesday, Oct 3
Santanu and I met again this morning. We agreed that with FTSM, it would be beneficial to
to support an array of multiple propagators and equilibrate them towards each individual phi psi
before actually updating free space.
Still however, the string is not correct. Eventually it becomes unstable.
Looking deeper at my equations...
Tuesday, Oct 2
Santanu is now saying we're back to a compiler error on HPCC again, seems it's tough to win here...
Fixed it in both places, I think this time for good.
Now moving back to the forces. I don't think I can straight up change the binding
of the nonbondedForceFunction data member. I'd have to somehow bind it to a Python function,
but the type of this function is from the template, and there are some other references to this
template data type. It could screw things up big time.
SWIG does however allow extension of a class; so I might be able to define a second evaluate() method
which instead of using the nonbondedForceFunction object of OneAtomPair, actually accepts
a nonbondedForceFunction and invokes an evaluate method with the forces.
I learned a valuable lesson with SWIG today! So apparently by default SWIG does not generate a
constructor for an abstract class. To get around this for an abstract class Foo:
%feature("notabstract") FooGood to know!
Except...for one problem....
It doesn't work. I get the same no constructor defined as before. Premature celebration. :(
Anyway, worked with Santanu some more today, we worked out the
bug in computePhiDihedral() and I reran FTSM. The string is definitely now tending in
the correct direction for alanine dipeptide. But it's not stable.
I think now might be a good time to put in the linear interpolation (as opposed to cubic
spline, which we talked about at Starbucks).
Put it there but I still get a string which collapses on itself if I do not fix the ends, which is not
what's supposed to happen.
Santanu is looking some more at the equations; I am helping Nava debug a ProtoMol problem..
Wrote her.
Now, I'm going to look at the force hooks again while Santanu checks the equations.
Monday, Oct 1
Working this morning on brainstorming a design for the forces...
Santanu came in and I'm going to discuss FTSM with him.
Okay, regarding the forces, within OneAtomPair there is an attribute
nonbondedForceFunction.
If I can somehow set this attribute to a force defined in Python.
I could inherit a new force from SystemForce.
Trying to think about how the control flow will work.
I think the ideal thing to do could initialize a force object based on the algorithm
(i.e. NonbondedCutoff, SimpleFull, PME, etc.) then I could set the nonbondedForceFunction
object to something that I create in Python.
Making an attempt at this...
Friday, Sept 28
Alrighty, today I should have some more time to try to figure out this obnoxious energy
problem with the Leapfrog Argon. Getting closer to solving it, I know the problem is with
the energy and not with the force...
My changes take a while to recompile since they're with a commonly used force, so in
between compilation I'm going to work on these epydoc comments.
Alright! Fixed the bug, the energies were not being cleared when the forces
were being cleared. Updated and e-mailed Dr. Izaguirre.
Now I'm spending this afternoon working on FTSM, the epydoc things are easier,
so I"ll save those for the weekend.
I have derivations for these derivatives, now I've just gotta figure out how to
use numpy to set up vectors, normalize them, etc. Must be something available.
Put in a computation of M and set up FTSM to run independent simulations.
Tried it and the results are still no good. But, we're probably further along than before
so I committed the code. It would be nice if I had some acceptable values of kappa and gamma.
Intermittendly adding more comments, would like to try out this epydoc soon, this
looks pretty neat.
Thursday, Sept 27
Alright, made my slides for tomorrow...
Now, starting to look at the epydoc which Dr. Izaguirre sent me.
This looks useful, the standard pydoc comments will generate little @ doxygen-like
comments, but if we generate good enough documentation with epydoc anyway
this shouldn't be too much of an issue.
I can start working on that, but right now I want to spend some time and
attempt a derivation of the equations for the M matrix.
Tonight I played with the leapfrog_Argon.
Figured out late that the total energy is not accurate. I know this
because I printed out the force vectors at the end of Leapfrog_Argon
and leapfrog_Argon and they are identical. So that means something
is not being set properly.
Anyway, 10:30, I'll look at it tomorrow.
TODO:
1. Fix the leapfrog_Argon bug
2. FTSM, matrix M (use Matlab)
3. Get the epydoc comments
Wednesday, Sept 26
Almost got this I/O working...the plots work from Leapfrog, but the outputs
are seg faulting for some reason.
Yep, that definitely is not right either. The object leapfrog consistently hangs
around 198 for the energy.
Alright I fixed the timestep, but the energy is still dropping.
Let me check its initial value...
Nope, very different (198-141).
Going to interdisperse this with some grading this afternoon..
Tuesday, Sept 25
More with MDL today...
Almost got things, just now a bug in the velocity reading.
Got it and recommitted.
Dr. Izaguirre pointed out a problem to me that propagator methods
do not run I/O.
Now that I moved I/O to be completely behind the scenes, having this loop
inside propagator methods is not working too well, because it requires an
explicit call to I/O in the method.
This is really not too good....
But the only other option is, don't pass numsteps to the propagator method
and don't use a loop, but that knocks out any ability to pre or post process
I'll wait for his reply on that, in the meantime I'm going to look into this FTSM
and making it work properly.
Fixed a few things, but I'm going to need some help from Santanu
on this M matrix, and also the last term of the free space update.
Like lambda(s,t), I have no idea what that is.
Looking at M....
M = sum(k=1, N)((1/m_k)*(dPhi/dx_k)*(dPsi/dx_k))
Thinking about this now, only a few atoms will have to be summed here:
the ones on which phi and psi depend.
Actually, only three atoms will contribute (the middle three).
Well now, but M is a function of x, so the equation is actually:
M(x)= sum(k=1, N)((1/m_k)*(dPhi(x)/dx_k)*(dPsi(x)/dx_k))
This doesn't make sense to me at all now...
Anyway, returning to io with the functions, I have to change the I/O
cache to not expect a propagator object now, which si fine...
Have to move Configuration and Integrator out of the Outputs.
Anyway, I have to shift gears for a minute, and get things together
for my class' exam tomorrow.
Monday, Sept 24
I hit my first hurdle with the FTSM code yesterday.
Basically, I build the topology for the system automatically when a user propagates
the system.
Now, we want to be able to access Phi and Psi before propagation. Which really, we
should be able to do. We should be able to get these the moment we set values in the
topology. Unfortunately, that's not how things work right now.
Way I see things, I have a few options. I could simply rebuild the topology every
time I set a related attribute, although that would be horribly inefficient it would
get the job done. I could also rebuild things everytime I access an attribute of the topology.
In a way, sort of a 'lazy' build. What I have to watch is, if I change boundary conditions,
I actually have to create a new topology object.
I think the latter makes more sense. One thing I must say I don't like about this
is that the user has to set boundary conditions before anything else.
Actually, I've got it. I can just initialize a default topology object, to Vacuum boundary conditions.
Then if they are reset, simply rebuild the topology again.
I have defaults for most things, so we should be safe here doing repeated builds.
K, after class...
I am having some strange problems with this computePhi() function. Somehow, no matter
which dihedral number I pass it always returns the exact same value for the angle.
Sunday, Sept 23
Looking at the FTSM code again this afternoon...
Thursday, Sept 20
Dr. Izaguirre sent me an e-mail that most of mdl was working.
I fixed a few more bugs last night, so hopefully now we are okay with
SWIG 1.3.31 and Python 2.5!
I am outlining the sets of goals that Dr. Izaguirre helped me outline
back in early August:
1. We need to build the hooks into energy, forces, and hessians
2. we need to efficiently support hessians (matrices) using ideally
numpy or existing library
3. we need to complete the FTSM implementation
4. we should compile a full set of examples on alanine dipeptide -
that allows us to use FTSM, but also other cool methods, and we can
compute interesting stuff on it. this could be very helpful for a
class. Examples of things to compute: free energy ramachadran plots
(easy, using matlab scripts already in use by our group); transition
rates (also easy, ditto); committors and tubes (new, needed by FTSM)
5. usability wise, we should:
a - have a Windows version
b - find an IDE
c - define a clean translation mech between protomol and mdl
d - test again interface to vmd/pmv
6. test performance
The clean translation between protomol and mdl has probably been
achieved (although could be improved) at this point.
I have to consider carefully if I want to produce all these tarballs yet,
or should I try to add some of the above first. Creating all these
tarballs I'm sure will be a day job.
As far as a "full set of examples" on alanine, basically we have
two right now, solvated and unsolvated.
The most urgent thing is probably to figure out the force hooks, although
it's also one of the most difficult.. the design of ProtoMol actually may
inhibit this slightly, I've gotta think about ways around...
The way evaluate() currently works for forces is it accepts:
topology, positions, velocities, forces and energies, then simply populates the
force vector.
Maybe this won't be too bad? The question then becomes, can I construct
a Python force (inheriting from Force), which in the background calls evaluate()?
Experimenting a little bit...my goal for now will be to add my own force to a MDL
ForceField, and get its evaluate() method called.
An immediate caveat I see is that evaluate() is called by the back end in evaluateSystemForces().
Therefore, it's not going to be as simple as adding a force that I've defined in Python
to a ForceField with an evaluate() that I wrote, and have it autmoatically called, because
the shared objects won't call the Python evaluate() but they'll attempt to call SystemForce::evaluate().
So what I may need to do is, similar to the propagators, be able to check the type
of force being defined, whether it is a Python force or an "optimized" SWIG-wrapped force
and handle them slightly differently, of course all behind the scenes.
Anyway, all they're doing is updating the force vector.
Which means I also need to ensure successful mapping to the back end force vector.
Fixed a bug, the forces were not being mapped to the back end properly...
Okay the difference here is this: with the propagators, I called initialize()
and run() from Python. So if I redefined them in Python we were okay, i.e.
if the propagator object was bound to a SWIG-wrapped propagator, the binding
of initialize() and run() would be to the wrapped versions, and if the propagator
object was a Python propagator, initialize() and run() would be bound to the Python
methods.
Here, I do not called evaluate() from Python. I call calculateForces() for Python
propagators, which loops through all forces and calls their evaluate() methods.
Making things worse, calculateForces() for SWIG-wrapped propagators is also
called from the wrapped propagator objects, so even this function if I redefined it
in Python would not be called for wrapped propagators such as BBK.
Wednesday, Sept 19
It's been a busy week for teh crouse, need to put together a study guide
and their first exam...
Monday, Sept 17
Yesterday and today, been working on compatibility with
the later version of Python/icc...
Python2.5 requires SWIG 1.3.31 or later...
There are good and bad to this, the good is this newer version
of SWIG simplifies template generation, the bad is that it
catches every memory leak. Although that could be considered good.
But still, I've got to remove every one of them before I can consider
MDL runnable. :p
Oh! There was a serious bug in the ./fullCopy.sh script. I was under
the assumption that >! overwrote a file in the output redirection. Turns out,
if the file exists it does nothing!!!
This fixed all of the memory leaks. W@W.
So now, I'm just dealing with a seg fault at the end of each run. Oh, and FTSM
seg faults at the end of the first step. This could have something to do with the seg fault at
the end though because of the individual independent propagations.
Got it, okay FTSM works and all other seg faults are gone.
Sent Dr. Izaguirre an e-mail after committing everything...
Now I have to work on the next assignment for class.
Friday, Sept 14
Committed the HPCC script finally. It works.
Doing some more tests with MDL today, making it robust...
Fixed several bugs in outputting data.
Thursday, Sept 13
Dr. Izaguirre e-mailed me this morning, I'm going to start looking
at this farm today.
Got the slides made, checking out the farm...
As expected, problems. Not letting me log in with the password I was provided.
Ah, okay.
Apparently you need to first create an account for their website,
and then create a second account to use the cluster. Was confusing.
Great, so I've got to wait another business day potentially...
Looking at hpcc.
Coded a new configure.in which statically links simtklapack, but I can't test
it because the queue is insane, so I just attached it in an e-mail to Dr.
Izaguirre.
Got the accoun ton the farm finally, now just need to figure out how to use it.
Wednesday, Sept 12
Class, slides online.. now the arc length reparameterization.
Fixed a few bugs, and now it updates the string repeatedly!
I am attempting to run an alanine system which I had been trying to get working
with MDL/FTSM using explicit solvent. I have removed the solvent
manually and run with ISM using CoulombDiElec, and after 20 ps of updates,
I get the green string in the figure below (the blue is the initial):
Going to see if there's someplace where I can verify this.
Tuesday, Sept 11
I am having weird and crazy deja vu this morning on Tues 9/11. Anyway..
Need slides slides slides for tomorrow...
AFS is down...again...
Finished slides.
Per a conversation today, Dr. Izaguirre is going to New York tomorrow.
My job for the rest of the day today is to try an FTSM script using
SCPISM (the Coulomb DiElec force).
Ok, I have a flat out simulation set up and running without explosions.
Now, looking at FTSM and setting things up again...
Monday, Sept 10
Ok. I have done some research and found that actually SWIG 1.3.29 and Python 2.5 are
noncompatible. There are some bugs.
However, SWIG 1.3.31 (the newest version) does work with Python 2.5 and I've been trying
to get linprotomol/mdl working with these guys. I have fixed Dr. Izaguirre's error,
but am getting new ones (the anomaly seems to be in which version of Vector3DBlock
is being used); also this swig/python combo is actually smart enough to detect for memory leaks,
and I had quite a few of them (esp. with the normal modes). So I'll probably have to get these
out before we release. Doesn't help that everything in Python is a reference. :p
Okay, good news - I have BBK working with 1.3.31 and 2.5, albeit with a seg fault at the end
that I"ll need to once again fix. :p
Man, this is incredible. So many things broke with these new versions.
K, we're getting close, fixed the memory leaks..
The only problems we have now are:
(1) Seg fault at the end
(2) Normal mode works with alanine but seg faults on WW
(3) HMC integrators have a memory leak at the end on a Vector3DBlock*, guessing
because they don't accept any forces.
(4) MDL integrators crash.
Fixed (4) but they now have a memory leak like (3).
Still lots of bugs as of almost 5, giving up for now will hit it again later.
Friday, Sept 7
Committed new binaries to sourceforge ProtoMol which statically
link the gcc library into the executable.
In the meantime, Dr. Thain sent me a link to a compile farm
available through NMI:
http://nmi.cs.wisc.edu/
So I'll take a look at that while testing a fresh checkout
of linprotomol/mdl from simtk.
Account is pending approval...
A fresh checkout just seemed to work out of svn,
so I'll have a look at Makefile.swig
Thursday, Sept 6
This morning, making up some more slides...
This afternoon, perfecting linprotomol by removing warnings in the compilation.
Once we are happy with things, I'll put up all those tarballs...
Removed warnings (some of them I had to suppress, because even though
SWIG would fail to wrap certain routines, they may still be needed in the back
end shared libraries).
Compiling, will recommit when done..
K, warnings are gone, seems to compile and copy over okay.
Looking at the machines that simtk lapack uses...
32-bit Linux AMD Generic
32-bit Linux AMD Opteron 2
32-bit Linux AMD Opteron 4
32-bit Linux AMD Opteron
32-bit Linux Intel Pentium 4-2
32-bit Linux Intel Pentium 4 -> manganese.helios.nd.edu (Development platform)
64-bit Linux AMD Opteron 2 -> lcls.cse.nd.edu, stones.cse.nd.edu
64-bit Linux AMD Opteron 4
64-bit Linux AMD Opteron -> steinbeck.helios.nd.edu
32-bit Mac G5-2
32-bit Mac G5-4
32-bit Mac G5 -> yahweh.cse.nd.edu
32-bit Windows AMD Generic -> doors.cse.nd.edu
32-bit Windows AMD Opteron 2
32-bit Windows AMD Opteron 4
32-bit Windows AMD Opteron
32-bit Windows Intel Generic
32-bit Windows Intel Pentium 4-2
32-bit Windows Intel Pentium 4
Now, I'm sure that I don't have all of these. Labelling above the ones
that I have access to. Sadly, it's a small percentage. :-(
ISS and biocomp fall into the above architectures too, although ISS is
an Intel Xeon.
Wednesday, Sept 5
Few more quickies this morning before class...
Done with class, back to MDL and wrapping this ProtoMol with g++.
Unfortunately my computer crashed yesterday, hopefully everything was saved. :P
Got NumPy on ISS, now fixing some errors with typenames that
icpc is fine with but g++ is not...
Ok, g++ works. Recommitting it, and verifying functionality with icpc...
Seems like we're okay with linear ProtoMol now. Good.
Now, interfacing things to MDL. Right now I have a ./fullCopy32.sh
and a ./fullCopy64.sh script. I think now we can remove all differentiation
within MDL between 32 and 64 bits. We can require people to install
Numpy, and they can also install Matplotlib. As far as I see, they're the only
two externals at the moment. Oh, right, SimTK Lapack. :p
Tuesday, Sept 4
All morning, writing the class slides and the next homework.
This afternoon, committed a linprotomol to svn.
Also working on the document Dr. Izaguirre sent me..
Done.
Trying to wrap linprotomol on iss, hoping the afs-dependent
errors will come up...
Okay, I'm running into my first hurdle, and it deals with a difference
between g++ and icpc it appears...
Value.h errors with g++. Doesn't like a forward declaration of Value
as being enough for a member routine to take Value& as an argument.
Lucky I checked with g++, it's spitting out a lot more linking errors than
icpc did.
Monday, Sept 3
Some more MDL stuff today...
YES!!! I am currently running a Normal Mode simulation on the ww domain
with MDL. Beautiful.
Now, that double free at the end is annoying me. Going to try to figure
that one out..
Done!
Okay, now I want to think about how to set up MDL on svn and improve
what we've got right now.
It'd be kind of nice to have my wrapped code in SVN. The only annoying
thing about that is if the user downloads MDL from the SVN repository, they'll get ProtoMol
as well.
Hmm, perhaps for now I could just keep my own copy of the wrapped ProtoMol
in my own cvs. This way I can release the source code as a tarball on the SVN site,
but still also have the ability for concurrent access on my own version.
So then the question becomes, which set of shared objects should I keep
in SVN?
Right now I have 32-bit and 64-bit Linux, plus Mac.
Looking at what SimTKLapack had...
Wow there are SO many.
Ooh, there will be an unforseen issue with the wrapped code and CVS.
How will I get new updates to Protomol and merge them in? Up till this
point I've had my wrapped code linked to SourceForge CVS so I could get
those updates.
Lol, hybrid SVN and CVS? I wonder if that's even possible. Don't see why not.
Trying to draw up on paper what this will look like.
Ok, so I'm envisioning the following:
On my development machine, I have a copy of, call it PyProtoMol. This stays
connected to SourceForge CVS so I can get all the recent updates.
Then, I'll set up this same directory in simtk SVN. On there is where I will commit
all updates, and retrieve updates across various architectures.
The only initial issue I see is that in general, to set up an SVN repository, you have to
first start a project, they'll give you a blank repository that you have
to check out, and you put your code there.
So the question will then become, will my CVS still work if I move my PyProtomol
to a new location? I think it will.
Then in PyProtoMol, I can have two scripts, one that builds (accepts a compiler
as an argument) and one that installs (accepts a path to MDL as an argument).
Fine.
Now the other issue is, the MDL on SVN. Currently I only have a specific set of
shared objects installed. It'd be nice to have all of them there I guess.
I wonder what they do with lapack? Like if I check out simtk lapack from the repository.]
Submitted linear ProtoMol for approval, awaiting...
In the meantime, trying to run in some more 'unfriendly' environments without AFS...
linprotomol approved tonight. Will start making it tomorrow morning.
Sunday, Sept 2
First thing today is to prepare lecture for Monday.
Then next up will be MDL. Since Normal Modes now seem to work with
Alanine, I want to
1. Test them on the WW
2. Fix that double free at the end
3. Commit new gcc to shared objects to svn, and maybe make some
downloadable gcc tarballs.
Lecture and homework done.
Friday, Aug 31
Last minute lecture changes this morning, then Normal Modes!
Ok, I found a descrepancy between MDL and ProtoMol. Unfortunately,
it deals with random numbers. Somehow, the first call to randomGaussian()
in NormModeInt::drift() returns different values between the two software packages. :(
I have checked and they do pass the same seed.
Somewhere, they are getting off track. Cannot find any correlation between them.
Ok, I've got the problem. It's that static int first that's causing trouble.
Got it ... now further, still some kind of bug...
Ok, I've found it I think.... thinking what I can do about it.
I am SO CLOSE now! There is still something small wrong, but the energies are
close (it still crashes with that Minimization Failed thingy, albeit after a much longer time).
Checking outputs...
Okay, there was one other bug in a velocity update, now we are extremely close, still
Alrighty!!! I've got matching output.
Now I've gotta get rid of these couts. Man there are so many, this could take some time
actually.
Thursday, Aug 30
I'm going to try to nail Normal Modes today once and for all...
While they compile, I'm testing out MDL with JC's installed Numpy
for 2.3.
Unfortunately, MDL is now completely broken with 2.3. It's going to
take me hours to fix this.
It seems that we have an urgent situation with MDL now... and that we're going
to have to do an overhaul again. :(
Users are just going to have to configure the language themselves.
I'll include the ProtoMol source. Sent Dr. Izaguirre a more detailed e-mail.
In the meantime, I'm going to keep working on Normal Modes
until I hear from him, because if I can get these working I can just pop
them into the Pm source that I include with the tarball.
Man, the forces in this simulation are different from the very beginning.
Wait, no, hang on... it looks like they're just being evaluated in different orders?
Hmm.
Yeah, in MDL for some reason Dihedral/Improper are evaluated before Angle/Bond.
So bonded forces actually match up.
LJ is off *slightly*, but the CoulombDiElec is way off (for atom 1's x component, we're
dealing with -0.04 vs. 1.92). So something is still seriously wrong here.
Trying a run with regular Coulomb...
AH, ok. I wasn't using exactly the same configuration parameters. NOW LJ lines up.
DiElec, however, still is way off. So that's the problem right there.
Ok, the switching function is working just fine now, it's the computed value
before the switch that is off.
Okay, I found the problem finally, the default value of D and s were off in the
constructor I was using for CoulombDiElec.
Alrighty!! NOW the forces match. w@w.
But my initial temperature is *way* off. (217-250). My initial energy is slightly
off, but I'm thinking it has something to do with this temperature.
Could be an issue with random numbers.
No, I've checked and the initial velocity values are all corresponding.
So it's got to be a difference in the initialize() section of the integrator. That's the only thing
I can think of at this point.
Okay! I have isolated things to a call to subspaceForce()!
mhQ appears different, let's see...
Okay, now fDof is 6 in MDL and 0 in ProtoMol... wow, why??
Fixed that problem... still however, that initial temperature is WAY off.
And so is the energy. Man, frustrating...
Ok, I lose the velocity correllation before NormModeInt gets called.
OHHH!! I wasn't removing angular momentum!! Just linear!!!
Whew. Ok step 0 is the same now. Man...
Unfortunately, minimization is still failing and by the end of the first
step output is radically different, so there's still more bugs. But, this is significant progress.
I now MUST switch to preparation for tomorrow's lecture.
Wednesday, Aug 29
Taught class and did office hours today; went further
with Normal Modes without much luck.
Dr. Izaguirre uncovered a serious MDL bug this afternoon
in the output ability, writing to files.
It seemed that all of them were broken because of
a back end change to the constructors! Ouch.
Fixed the 64-bit, 32 and Mac still need fixing.
This is why we need a fully configurable version.
Fixing 3 versions is a pain.
Tuesday, Aug 28
Running normal mode tests this morning, updating my PM baseline...
The file I am running is alan.nm.conf. Output I get is here:
Step : 0, Time : 0.000 [ps], TE : 5.0750 [kcal/mol], T : 217.0237 [K], V : 168.69 [AA^3]
Step : 100, Time : 0.400 [ps], TE : 12.6803 [kcal/mol], T : 365.0835 [K], V : 171.70 [AA^3]
Step : 200, Time : 0.800 [ps], TE : 11.5982 [kcal/mol], T : 299.7578 [K], V : 152.91 [AA^3]
Step : 300, Time : 1.200 [ps], TE : 9.3692 [kcal/mol], T : 243.7300 [K], V : 176.44 [AA^3]
Step : 400, Time : 1.600 [ps], TE : 15.0602 [kcal/mol], T : 414.7597 [K], V : 164.38 [AA^3]
Step : 500, Time : 2.000 [ps], TE : 11.7877 [kcal/mol], T : 253.9186 [K], V : 150.31 [AA^3]
Step : 600, Time : 2.400 [ps], TE : 8.2995 [kcal/mol], T : 240.8821 [K], V : 164.16 [AA^3]
Step : 700, Time : 2.800 [ps], TE : 13.9773 [kcal/mol], T : 325.3749 [K], V : 170.11 [AA^3]
Step : 800, Time : 3.200 [ps], TE : 11.9056 [kcal/mol], T : 394.3018 [K], V : 176.52 [AA^3]
Step : 900, Time : 3.600 [ps], TE : 11.3224 [kcal/mol], T : 371.0887 [K], V : 207.87 [AA^3]
Step : 1000, Time : 4.000 [ps], TE : 9.4503 [kcal/mol], T : 351.1716 [K], V : 195.79 [AA^3]
Step : 1100, Time : 4.400 [ps], TE : 9.8656 [kcal/mol], T : 312.4039 [K], V : 196.81 [AA^3]
Step : 1200, Time : 4.800 [ps], TE : 9.6716 [kcal/mol], T : 252.3893 [K], V : 171.44 [AA^3]
Step : 1250, Time : 5.000 [ps], TE : 14.7610 [kcal/mol], T : 477.4991 [K], V : 150.50 [AA^3]
Hint: [NormModeInt::~NormModeInt] Minimizer iterations = 1, average = 1.2032, force calcs = 3, average = 3.4064
Timing : wall: 1.23864[s] real, 0.85487[s] user, 0.08099[s] sys, run: 0.89194[s], integration: 0.69901[s], forces: 0.55692[s].
MDL output (blah):
Step : 0, Time : 0.000 [fs], TE : 11.7903 [kcal/mol], T : 250.8767 [K], V : 168.69 [AA^3]
python: ./framework/base/Array_Fastest.h:453: ProtoMol::RefArray<T, <expression>>
&ProtoMol::Array<T, N>::operator[](unsigned int) [with T = std::pair<ProtoMol::CubicCellLocation, int>, N = 3U]:
Assertion `Index<m_NDimensions[0]' failed.
Bugs bugs bugs.
This DiElectric coulombic force is still not right. I get crazy values when I add that to
the forces to evaluate... sigh, frustrating.
Okay, it's definitely the switching function. But I thought I debugged this!!! Grr.
I got the switching function problem fixed, at least now velocities/forces are acceptable.
Alas, still I'm getting the above crash:
Step : 0, Time : 0.000 [fs], TE : 1.4405 [kcal/mol], T : 250.8767 [K], V : 168.69 [AA^3]
python: ./framework/base/Array_Fastest.h:453: ProtoMol::RefArray<T, <expression>>
&ProtoMol::Array<T, N>::operator[](unsigned int) [with T = std::pair<ProtoMol::CubicCellLocation, int>, N = 3U]:
Assertion `Index<m_NDimensions[0]' failed.
So something is still wrong. About step 22, I'm getting nan positions.
Fixed another bug, but now I get at step 25:
Fatal Error: [NormModeMin::updtConstrainedModes] Minimization failed. Rediagonalization maybe required for
the current conformation. Aborting.
Still stuff going on.
Anyway, Dr. Izaguirre is wanting a configurable version of MDL. It'll take a while,
I'm going to map out a plan of attack tonight and tomorrow after my class.
Monday, Aug 27
Figured out that it's the switching that's causing these nan and inf LennardJones forces
with Normal Modes...
Wow, I found the problem. I was setting the cutoff too late in the forces, after other
parameters were being set which depended on the cutoff. Therefore, when I setup the forces
the cutoff always needs to be before the extra paramters.
My fingers are crossed now...
ALRIGHT!!! SQUASHED IT!
Now, gotta clean up my back end stuff...
Ok, now it runs, so that was a huge bug fix. But I do still get exploding energies, although
this time I have a suspicion as to why. ;-)
Yep. Better now, but still eventually the energy gets too high and I get this message:
Hint: [NormModeMin::updtConstrainedModes] Reset CG, PE fail. rsCG= 1
Hint: [NormModeMin::updtConstrainedModes] Reset CG, PE fail. rsCG= 2
Hint: [NormModeMin::updtConstrainedModes] Reset CG, PE fail. rsCG= 3
Hint: [NormModeMin::updtConstrainedModes] Reset CG, PE fail. rsCG= 4
Hint: [NormModeMin::updtConstrainedModes] Reset CG, PE fail. rsCG= 5
Fatal Error: [NormModeMin::updtConstrainedModes] Minimization failed. Rediagonalization maybe required for the current conformation. Aborting.
So still more problems. :(
It may be that my input files are out-of-date as well. E-mailed Chris...
From what I can see in these simulations, it is without a doubt the kinetic
energy that is skyrocketing.
Which means, probably the velocities are getting out of whack.
Okay, I am noticing one suspicious thing. .. my coulombic energy is zero.
I don't think that should be so. Perhaps not accounting for that value is throwing
everybody else off.
No, I was wrong about this... I was indexing the scalarstructure wrongly. We're fine
with Coulombic.
Santanu just sent me a bunch of input files that he knows work with NM. I'm going
to try those...it's possible these WW ones are dated.
No, I've tried these, gives me the above error right away.
GRRRRRRRRR!
The velocities in the early going of this run are all nan. So definitely, I am killing these
velocities somehow. But where, or where have my velocities gone?
Okay, I wasn't reading the eigenvector file properly.
It's better now, but still crashes eventually with the above errors. :(