Dear All, I am running Langevin dynamics on a coarse-grained protein model using Charmm (version c35b5, single CPU). I've calculated some time-dependent correlation functions and observe periodic correlations at specific time-intervals that are completely unnatural and should not occur. I've attached a graph showing one such correlation function, with the periodic correlations indicated by red arrows. The correlation function is the auto-correlation of a protein domain vector over time.
Do these unnatural periodic correlations in the graph indicate I have exceeded the period of the random number generator used in the random force generator of Langevin dynamics? I imagine exceeding its period would result in correlated forces over long times and result in this observation.
Is there any way to get around this issue? I would like to run long simulations without such artifacts.
Note, that I am using the NEWRNG random number generator in the simulations and I apply Langevin dynamics to 171 interaction sites with a damping coefficient (fbeta) of 5. The periodic repetition in the correlation function occurs every 890,000 ps of simulation time. Therefore, I estimate that at this point in the simulation NEWRNG has been called at least 2.28 x 10^9 times (= 890,000 ps * 5 ps^-1 * 171 interaction sites * 3 force components).
Any comments / suggestions about the origin of this issue and how to get around it would be appreciated. Thanks, Ed
Here is output demonstrating I use NEWRNG: ~/projects/>$c35b5_dhdwp < temp.inp 1 Chemistry at HARvard Macromolecular Mechanics (CHARMM) - Developmental Version 35b5 August 15, 2010 Copyright(c) 1984-2001 President and Fellows of Harvard College All Rights Reserved Current operating system: Linux- Created on 3/20/11 at 9:43:59 by user: test
Maximum number of ATOMS: 360720, and RESidues: 120240 Current HEAP size: 10240000, and STACK size: 10000000
RDTITL> * RDTITL> No title read.
***** LEVEL 1 WARNING FROM <RDTITL> ***** ***** Title expected. ****************************************** BOMLEV ( 0) IS NOT REACHED. WRNLEV IS 5
Here is my CHARMM script to run Langevin dynamics:
!! Set up temperature run set timestp 0.015 set steps 320000000 set temp 312 set tstruc 312 set psf ../../setup/rnc.psf set param ../../setup/rnc.prm set top ../../setup/rnc.top set fbsolu 5.0 set strt setup/1_90000_prod.cor set OUT ddfln
! read parameter and topology files open unit 10 read form name @top read rtf unit 10 card close unit 10 open unit 10 read form name @param read param unit 10 card flex close unit 10
! Get psf read psf file name @psf SHAKE BOND PARA
! fix some atoms, and restrain linker to origin cons fix sele segid RIB end MMFP GEO sphere harm - xref 0 yref 0 zref 0 - force 50 droff 0.0 select resid @resll .and. segid A .and. type A end END
prnlev @prlev ! Get the coordinates open read unit 3 file name @strt read coor unit 3 file close unit 3
! speeds up the non-bonded list generation nbond bycc eten on
mini abnr nstep 1000
! Set friction forces scalar fbeta set 0.05 sele segid A end SCAL FBETA SHOW sele segid A end
The period of the new random number generator is about 2^121. I'm not sure you actually get the new generator in c35b5 unless you specifically ask for it: random clcg
With Langevin dynamics random numbers are generated in each integration step, not every 1/beta ps.
Lennart Nilsson Karolinska Institutet Stockholm, Sweden
Hi Prof. Nilsson, Thanks for the info. With forces included at each timestep it should only increase the number of calls to the random number generator from 10^9 to 10^12 or so. I was unaware that I had to specifically invoke it in the CHARMM script to get the NEWRNG generator:
random clcg ! Set friction forces scalar fbeta set @fbsolu sele segid A end SCAL FBETA SHOW sele segid A end open unit 21 write unform name output/@OUT_@traj.dcd open unit 22 write form name output/@OUT_@traj.mon dyna leap langevin strt nstep @steps timestep @timestp - iprfrq 0 ieqfrq 0 ntrfrq 0 - iunwri -1 iuncrd 21 iunvel -1 kunit 22 - ihtfrq 0 teminc 0 nprint 1500 nsavc 1500 isvfrq 0 nsavv 0 ihbfrq 0 - inbfrq -1 imgfrq 0 - ilbfrq 0 - ! and step frequency for writing restart file rbuffer 0.0 tbath @temp tstruc @tstruc - firstt @temp finalt @temp - iasors 0 iasvel 1 iscvel 0 ichecw 0 twindh 0.0 twindl 0.0 iseed @rand - echeck 1000
I'm running a test trajectory now with this correction to see if it avoids these problems.
Is it correct that I still need to invoke 'NEWRNG' in pref.dat? Are there any other keywords I need to include to access 'rand clcg'? Thank you very much for your help, Ed
the ##IF NEWRNG directive is no longer present in c35b? versions, so it has no effect
the CLCG random numbers (NEWRNG) are the default and are used for Langevin, Monte Carlo, etc.
the old IMSL random number code is used only for ?RAND
the command "RANDOM CLCG" sets logical QOLDRNG = .false., but QOLDRNG is never used elsewhere in the source code.
A weakness in this was that the seed values were not stored in the restart file, so DYNA LANG RESTART *requires* a new seed with c35b? versions of CHARMM. For the c36 versions that was fixed, and a lot of changes in random number behavior were made as well, but that code is not yet publicly available.
The first executable statements in function RANDOM in c35b3/source/util/clcg.src are: if (qoldrng) then random = oldrandom(g) return endif
and source/charmm/iniall.src has:
QOLDRNG=.TRUE. if (.not.qoldrng) then call CLCGInit(ISEED) ISEED=1 endif
I think the "default", as far as the new random number generator in c35 is concerned, is that the code is compiled by default; it does however not seem to be used by default.
I don't have c35b5 online, and it may be different. Ed is wise to check.
Lennart Nilsson Karolinska Institutet Stockholm, Sweden
For completeness. The period of the old random number generator is 2^31 - 2. It generates all integer numbers between 1 and this number exactly once in a pseudorandom order before it starts over again.
Lennart Nilsson Karolinska Institutet Stockholm, Sweden
c34bN The RNG keyword to install.com adds NEWRNG keyword to pref.dat, using the new CLCG RNG for all simulation methods; the seed values are not stored in the restart file
c35bN NEWRNG keyword removed, the CLCG code is included by default, but must be enabled via RANDOM CLCG command; seeds not stored in restart file
c36aN The CLCG RNG is the default, and the seed values may be stored in the restart file
