Previous Thread
Next Thread
Print Thread
Page 1 of 2 1 2
Joined: Aug 2004
Posts: 192
E
Edo Offline OP
Forum Member
OP Offline
Forum Member
E
Joined: Aug 2004
Posts: 192
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


CHARMM>

CHARMM> PREF

KEYWORD LIST: Number of keys: 109.
GNU UNIX SCALAR I4BINARY INTEGER8 XXLARGE EXPAND PUTFCM ACE
ADUMB AFM ASPENER ASPMEMB BANBA BLOCK CFF CHEQ CMAP DIMB DMCONS
DOCK DYNVV2 EMAP ESTATS FASTEW FLEXPARM FLUCQ FMA FOURD FSSHK
GBBLCK GBFIXAT GBIM GBINLINE GBMV GBMVFAST HDGB HQBM GBSW GBSWIT
GCMC GENETIC GENBORN GOMODEL GRID GSBP HMCM IMCUBES LATTICE LDM
LDLAN LDMGEN LONEPAIR LMC LRST LRVDW MC MEHMC MMFF MOLVIB
MULTCAN NBIPS NEWTIMER NEWRNG OLDDYN OPLS OVERLAP PATHINT PBEQ
PBOUND PERT PHMD POLAR PSSP PM1 PMEPLSMA PNOE PRIMSH RCFFT
RDFSOL REPLICA RGYCONS RPATH RXNCOR SASAE SCPISM SGLD SHAPES
SHELL SOFTVDW TNPACK TPS TRAVEL TSM TSALLIS FITCHG WCA CHEMPERT
PROTO SMD COMP2 FACTS LOOKUP RMD RXNCONS TAMD QUANTUM RUSH
NOGRAPHICS


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

! equilibrate
dyna leap langevin strt nstep 16000000 timestep @timestp -
iprfrq 0 ieqfrq 0 ntrfrq 0 -
iunwri -1 iuncrd -1 iunvel -1 kunit -1 -
ihtfrq 0 teminc 0 nprint 4000 nsavc 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 100000

! Set friction forces
scalar fbeta set @fbsolu sele segid A end
SCAL FBETA SHOW sele segid A end

! production run
calc rand @rand+7
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

Attached Images
corr_funct_vs_time.jpg (566.64 KB, 305 downloads)
Last edited by Edo; 03/20/11 02:54 PM.
Joined: Sep 2003
Posts: 4,863
Likes: 10
Forum Member
Online Content
Forum Member
Joined: Sep 2003
Posts: 4,863
Likes: 10
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
Joined: Aug 2004
Posts: 192
E
Edo Offline OP
Forum Member
OP Offline
Forum Member
E
Joined: Aug 2004
Posts: 192
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

Joined: Sep 2003
Posts: 8,629
Likes: 24
rmv Online Content
Forum Member
Online Content
Forum Member
Joined: Sep 2003
Posts: 8,629
Likes: 24
According to the source code--
  • 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.


Rick Venable
computational chemist

Joined: Sep 2003
Posts: 4,863
Likes: 10
Forum Member
Online Content
Forum Member
Joined: Sep 2003
Posts: 4,863
Likes: 10
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
Joined: Aug 2004
Posts: 192
E
Edo Offline OP
Forum Member
OP Offline
Forum Member
E
Joined: Aug 2004
Posts: 192
Yes Prof. Nilsson, c35b5 has the same corresponding source code (excerpted here):
if (qoldrng) then
random = oldrandom(g)
return
endif

QOLDRNG=.TRUE. ! .FALSE. !yw 05-Aug-2008
if (.not.qoldrng) then
call CLCGInit(ISEED)
ISEED=1
endif

I'm hoping NEWRNG (CLCG) isn't used by default in c35b5, otherwise I'll have to find a different source to my original error.

Joined: Sep 2003
Posts: 8,629
Likes: 24
rmv Online Content
Forum Member
Online Content
Forum Member
Joined: Sep 2003
Posts: 8,629
Likes: 24
I should probably remember to use the '-i' flag of grep ...


Rick Venable
computational chemist

Joined: Sep 2003
Posts: 4,863
Likes: 10
Forum Member
Online Content
Forum Member
Joined: Sep 2003
Posts: 4,863
Likes: 10
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
Joined: Sep 2003
Posts: 8,629
Likes: 24
rmv Online Content
Forum Member
Online Content
Forum Member
Joined: Sep 2003
Posts: 8,629
Likes: 24
If I have things correct now, the situation is--
  • 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


Rick Venable
computational chemist

Joined: Sep 2003
Posts: 4,863
Likes: 10
Forum Member
Online Content
Forum Member
Joined: Sep 2003
Posts: 4,863
Likes: 10
Shouldn't the third bullet "may be" read "are"?


Lennart Nilsson
Karolinska Institutet
Stockholm, Sweden
Page 1 of 2 1 2

Moderated by  lennart, rmv 

Link Copied to Clipboard
Powered by UBB.threads™ PHP Forum Software 7.7.5
(Release build 20201027)
Responsive Width:

PHP: 7.3.31-1~deb10u2 Page Time: 0.010s Queries: 36 (0.006s) Memory: 0.7910 MB (Peak: 0.9094 MB) Data Comp: Off Server Time: 2023-01-27 08:09:21 UTC
Valid HTML 5 and Valid CSS