Skip to content

Commit 7b31d2d

Browse files
ceriottmceriottm
committed
Added the Fortran example to the master branch
Trying to pull in as much material as possible from the old website.
1 parent da95a51 commit 7b31d2d

File tree

9 files changed

+753
-0
lines changed

9 files changed

+753
-0
lines changed

‎f90-example/GLE-A

+6
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,6 @@
1+
4
2+
3.0770109810000013 2.451344375000001 -0.004620728182000001 0.049762613900000016 0.3364690971000002
3+
2.451344449 2.5484052860000013 0. 0. 0.
4+
-0.004598788959000001 0. 1.7311834200000007 0. 0.
5+
-0.04976420041000001 0. 0. 0.9243109136000003 0.
6+
0.33646911120000006 0. 0. 0. 0.2728914538000001

‎f90-example/Makefile

+20
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,20 @@
1+
.PHONY: all clean
2+
FC := gfortran -x f95-cpp-input
3+
LD := gfortran
4+
FFLAGS := -O3 -mtune=native -DUSELIBS
5+
OBJS := $(patsubst %.f90,%.o,$(wildcard *.f90))
6+
7+
all: harmonic.x
8+
9+
harmonic.x: $(OBJS)
10+
$(LD) $(FFLAGS) -o harmonic.x $(OBJS) -lblas -llapack
11+
12+
md-nose.o: md-tools.o
13+
md-gle.o: md-tools.o
14+
harmonic.o: md-tools.o md-nose.o md-gle.o
15+
16+
%.o: %.f90
17+
$(FC) $(FFLAGS) -c -o $@ $<
18+
19+
clean:
20+
rm *~ *.o *.mod harmonic.x

‎f90-example/README

+38
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,38 @@
1+
* * *
2+
Generalised Langevin Equation - A toy code
3+
* * *
4+
5+
An awkwardly simple code to test Generalised Langevin
6+
dynamics on a multidimensional harmonic oscillator.
7+
Licensed under GPLv3 [www.gnu.org], source can be used
8+
as a stub to implement GLE thermostat in existing codes.
9+
For comments and requests, write to
10+
michele dot ceriotti at gmail dot com
11+
12+
13+
* Program files:
14+
harmonic.f90 main, input parsing and md loop
15+
md-tools.f90 a couple of utility functions
16+
md-nose.f90 simple implementation of NH chains
17+
md-gle.f90 GLE code, mostly self-contained
18+
19+
* Compilation:
20+
A Makefile is provided, even if the code is so simple
21+
compilation should be trivial. Adjust environment
22+
variables in the Makefile's header, and just "make".
23+
24+
* Use:
25+
Input is detailed in the source. Beside the actual input,
26+
one is required a 'freq' file containing the square root
27+
of the Hessian of the harmonic potential, and, if GLE is
28+
to be used, a GLE-A file containing the drift matrix,
29+
in the format
30+
ns
31+
app ap_1 ... ap_ns
32+
bap_1 A11 ... A1ns
33+
.
34+
bap_ns Ans1 ... Ansns
35+
If non-FDT dynamics is to be performed, a similar GLE-C
36+
is required, containing the covariance matrix.
37+
38+

‎f90-example/freq

+2
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,2 @@
1+
0.1 0.0
2+
0.0 0.01

‎f90-example/harmonic.f90

+183
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,183 @@
1+
! *********************************************************
2+
! * An awkwardly simple code for the dynamics of an *
3+
! * an n-dimensional harmonic oscillator, to demonstrate *
4+
! * the use of colored-noise Langevin equation in *
5+
! * thermostatting. *
6+
! * *
7+
! * Code has been kept as modular as possible, and it *
8+
! * should be relatively straightforward to adapt it for *
9+
! * other serial MD codes. In parallel codes using domain *
10+
! * decomposition, the additional degrees of freedom *
11+
! * should "follow" the corresponding atoms. *
12+
! * *
13+
! * Feel free to copy, modify, lend, borrow, nuke this *
14+
! * code, which is licensed under GPLv3 [www.gnu.org] *
15+
! * e-mail me at michele dot ceriotti at gmail dot com *
16+
! *********************************************************
17+
18+
program harmonic
19+
use md_tools
20+
use md_nose
21+
use md_gle
22+
23+
implicit none
24+
25+
! * options to be given in input file *
26+
! seed initial seed for PRNG
27+
! temp target temperature
28+
! dt timestep
29+
! nstep number of steps to be performed
30+
! stride output every stride frames
31+
! ndof dimensionality of the oscillator
32+
! wfile name of the file to read the hessian
33+
! decomposition from (see below)
34+
! traj logical, whether to output p,q trajectory
35+
! * thermostatting options *
36+
! wnw "optimal" frequency for langevin (0. 0-> WN off)
37+
! glew range shift freq. for GLE (0. 0-> GLE off)
38+
! nchains number of NH chains to be used (0 -> NHC off)
39+
! nhmts timestep subdivision for NHC integration
40+
! nhw "optimal" frequency for NHC (Q=kt/nhw**2)
41+
! * a note on thermostats: all the thermostats can be used at once,
42+
! even if that makes not much sense. again, this is just a small
43+
! test code
44+
! * a note on units: mass is set to one, and k_B as well, so that
45+
! at equilibrium <p^2>=omega^2 <q^2>=temp
46+
47+
real*8 dum, dt, temp, nhw, wnw, glew
48+
integer argc, seed, nstep, stride, ndof, nchains, nhmts, tstride
49+
character *256 fname, wfile, prefix
50+
namelist /inp/ seed, wfile, dt, temp, nstep, stride, tstride, ndof, nchains, nhmts, nhw, wnw, glew
51+
52+
real*8, allocatable :: q(:), p(:), f(:), wm(:,:)
53+
real*8 v, k, h ! yes, it's all we need!
54+
real*8 dt2
55+
#ifdef USELIBS
56+
integer, allocatable :: ipiv(:) ! needed to invert matrix
57+
#endif
58+
integer irnd, istep, i, j
59+
60+
! reads command line
61+
argc = iargc()
62+
if ( argc .ne. 1 ) then
63+
write(6,*) '* Call me as: harmonic <input>'
64+
stop
65+
endif
66+
67+
! reads input file
68+
tstride=0
69+
call getarg (1,fname)
70+
open(101,file=fname)
71+
read(101,inp)
72+
close (unit=101)
73+
irnd=-seed
74+
dt2=dt*0.5
75+
kt=temp
76+
ndim=ndof
77+
allocate(q(ndof))
78+
allocate(p(ndof))
79+
allocate(f(ndof))
80+
81+
allocate(wtw(ndof,ndof))
82+
allocate(wm(ndof,ndof))
83+
#ifdef USELIBS
84+
allocate(ipiv(ndof))
85+
#endif
86+
87+
! init random seed
88+
dum=ran2(irnd)
89+
90+
! reads square root of the hessian (v(q)=1/2 q^T W^T W q)
91+
open(101,file=wfile)
92+
read(101,*) wm
93+
close (unit=101)
94+
wtw=matmul(transpose(wm),wm)
95+
96+
! init momenta
97+
do i=1,ndof
98+
p(i)=rang(irnd)*sqrt(temp)
99+
enddo
100+
! init coordinates (requires inverting wm ONLY WORKS IF WM IS SQRT(HESSIAN), NOT IF DONE WITH CHOLESKY!!!)
101+
#ifdef USELIBS
102+
do i=1,ndof
103+
q(i)=rang(irnd)*sqrt(temp)
104+
enddo
105+
call dgesv(ndof, 1, wm, ndof, ipiv, q, ndof, i)
106+
#else
107+
! define your own inversion if you wish.
108+
! I am lazy and just init to zero.
109+
q=0.
110+
#endif
111+
call force(q,f)
112+
113+
!initializes thermostats
114+
if(nchains .gt. 0) call nhc_init(nchains, nhw, irnd)
115+
if(glew .gt. 0.d0) call gle_init(dt2,glew,irnd)
116+
if(wnw .gt. 0.d0) call wn_init(dt2,wnw)
117+
118+
! opens file for output
119+
if (tstride>0) open(102,file='traj-p.out')
120+
if (tstride>0) open(103,file='traj-q.out')
121+
open(104,file='statis.out')
122+
123+
! we are already at the main dynamics loop!
124+
do istep=1,nstep
125+
126+
!calls the active thermostats
127+
if(nchains .gt. 0) call nhc_step(p,dt2,nhmts)
128+
if (wnw .gt. 0.d0 .or. glew .gt. 0.d0) then
129+
call kin(p, k)
130+
langham=langham+k
131+
if (wnw .gt. 0.d0) call wn_step(p,irnd)
132+
if (glew .gt. 0.d0) call gle_step(p,irnd)
133+
call kin(p, k)
134+
langham=langham-k
135+
endif
136+
137+
! hamiltonian step for dynamics
138+
p=p+f*dt2
139+
q=q+p*dt
140+
call force(q,f)
141+
p=p+f*dt2
142+
143+
! thermostats, second bit
144+
if (wnw .gt. 0.d0 .or. glew .gt. 0.d0) then
145+
call kin(p, k)
146+
langham=langham+k
147+
if (wnw .gt. 0.d0) call wn_step(p,irnd)
148+
if (glew .gt. 0.d0) call gle_step(p,irnd)
149+
call kin(p, k)
150+
langham=langham-k
151+
endif
152+
if(nchains .gt. 0) call nhc_step(p,dt2,nhmts)
153+
154+
! computes properties & outputs
155+
if (mod(istep,stride).eq.0) then
156+
call pot(q, v)
157+
call kin(p, k)
158+
h=k+v
159+
if (nchains.gt.0) then
160+
call nhc_cons()
161+
h=h+nhcham
162+
end if
163+
if (wnw .gt. 0.d0 .or. glew .gt. 0.d0) then
164+
h=h+langham
165+
end if
166+
if (tstride.gt.0 .and. mod(istep,tstride).eq.0) then
167+
write(102,'(1e17.8 )', advance='NO') istep*dt
168+
write(103,'(1e17.8 )', advance='NO') istep*dt
169+
do i=1,ndof
170+
write(102,'( 1e17.8 )', advance='NO') p(i)
171+
write(103,'( 1e17.8 )', advance='NO') q(i)
172+
enddo
173+
write(102,*) ""
174+
write(103,*) ""
175+
endif
176+
write(104,'(4e17.8)') istep*dt, v/ndof, k/ndof, h/ndof
177+
endif
178+
enddo
179+
180+
if (tstride>0) close(102)
181+
if (tstride>0) close(103)
182+
close(104)
183+
end program harmonic

‎f90-example/input

+15
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,15 @@
1+
&inp
2+
seed = 12345
3+
temp = 1
4+
dt = 5e-2
5+
nstep = 2000000
6+
stride = 10
7+
wnw = 0
8+
glew = 1.0
9+
nhw = 0.0
10+
nchains = 0
11+
nhmts = 6
12+
wfile = 'freq'
13+
ndof = 2
14+
&end
15+

0 commit comments

Comments
 (0)