-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathQfep.f90
209 lines (154 loc) · 8.41 KB
/
Qfep.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
PROGRAM Qfep
! #DES: Program to produce and write the FEP/US calculation data in Qfep format. (SIC throughout)
USE Util, ONLY : Startup, Cleanup
USE Data, ONLY : ComputeDerivedData
USE Log, ONLY : logUnit
IMPLICIT NONE
INTEGER, PARAMETER :: outUnit = 64
CALL Startup(readCoords=.FALSE.)
CALL ComputeDerivedData(logUnit,doTiming=.FALSE.,readCoords=.FALSE.,doFEPUS=.TRUE.)
CALL QFepAnalysis()
CALL CleanUp()
CONTAINS
!Run the 3-stage analysis of the simulation data as in qfep
SUBROUTINE QfepAnalysis()
! #DES: Master routine: create a Qfep log file, do the analyses and write the output.
USE FileIO, ONLY : OpenFile, CloseFile
USE Log, ONLY : logUnit
IMPLICIT NONE
WRITE(logUnit,'(A)') "Performing qfep analysis - output to QFEPstyle.log"
WRITE(logUnit,*)
CALL OpenFile(outUnit,"QFEPstyle.log","write")
WRITE(outUnit,'(A)') "Qfep-style output generated by Fepcat (NB: small numerical differences with Qfep output may be observed)"
WRITE(outUnit,*) ""
CALL WriteParameters()
CALL AnalyzeDynamics_Q()
CALL AnalyzeFEP_Q()
CALL AnalyzeFepUs_Q()
CALL CloseFile(outUnit)
END SUBROUTINE QfepAnalysis
!*
SUBROUTINE WriteParameters()
USE Input, ONLY : nFepSteps, stateA, stateB, nStates, nSkip, nBins, minPop, rcCoeffA, rcCoeffB, alpha, beta
IMPLICIT NONE
WRITE(outUnit,'(A27)') "--> Number of energy files:"
WRITE(outUNit,'(A20,I4)') "# Number of files = ", nFepSteps
WRITE(outUnit,'(A55)') "--> No. of states, no. of predefined off-diag elements:"
WRITE(outUnit,'(A21,I4)') "# Number of states = ", nStates
WRITE(outUnit,'(A36,I4)') "# Number of off-diagonal elements = ", 0
WRITE(outUnit,'(A54)') "--> Give kT & no, of pts to skip & calculation mode:"
WRITE(outUnit,'(A7,F7.3)') "# kT = ", 1.0d0 / beta
WRITE(outUnit,'(A34,I6)') "# Number of data points to skip = ", nSkip
WRITE(outUnit,'(A44,L1)') "# Only QQ interactions will be considered = ", .FALSE.
WRITE(outUnit,'(A28)') "--> Give number of gap-bins:"
WRITE(outUnit,'(A23,I4)') "# Number of gap-bins = ", nBins
WRITE(outUnit,'(A27)') "--> Give minimum # pts/bin:"
WRITE(outUnit,'(A37,I4)') "# Minimum number of points per bin = ", minPop
WRITE(outUnit,'(A27)') "--> Give alpha for state 2:"
WRITE(outUnit,'(A22,F7.2)') "# Alpha for state 2 = ", alpha(stateB) - alpha(stateA)
WRITE(outUnit,'(A16)') "--> Hij scaling:"
WRITE(outUNit,'(A25,F7.2)') "# Scale factor for Hij = ", 1.0d0
WRITE(outUnit,'(A57)') "--> linear combination of states defining reaction coord:"
WRITE(outUnit,'(A37,2F7.2)') "# Linear combination co-efficients = ", rcCoeffA, rcCoeffB
WRITE(outUnit,*); WRITE(outUnit,*)
END SUBROUTINE WriteParameters
!*
SUBROUTINE AnalyzeDynamics_Q()
! #DES: Reproduce the output simulation data table of Q (same data, slightly different format)
USE Input, ONLY : energyNames, stateEnergy, coeffs, mask, CreateFileNames, fileBase, nFepSteps
USE StatisticalFunctions, ONLY : mean
IMPLICIT NONE
INTEGER :: name, fepstep, state, type
CHARACTER(500) :: fileNames(nFepSteps)
CALL CreateFileNames('en',fileBase,fileNames)
WRITE(outUnit,'(A)') "# Part 0: Average energies for all states in all files"
WRITE(outUnit,'(A9,1X,A5,1X,A8,1X,A9)',ADVANCE='NO') "FEP Index", "State", "Points", "Coeff"
DO name = 1, SIZE(energyNames) !nEnergyTypes
WRITE(outUnit,'(A11)',ADVANCE='NO') energyNames(name)
ENDDO
WRITE(outUnit,*)
DO fepstep = 1, SIZE(stateEnergy,2) !nFepSteps
DO state = 1, SIZE(stateEnergy,3) !nStates
WRITE(outUnit,'(A,1X,I5,1X,I8,1X,F9.6)',ADVANCE='NO') TRIM(ADJUSTL(fileNames(fepstep))), state, COUNT(mask(fepstep,:)), coeffs(1,fepstep,state)
DO type = 1, SIZE(energyNames)
WRITE(outUnit,'(F11.2)',ADVANCE='NO') mean( stateEnergy(:,fepstep,state,type),mask=mask(fepstep,:) )
ENDDO
WRITE(outUnit,*)
ENDDO
WRITE(outUnit,*)
ENDDO
ENDSUBROUTINE AnalyzeDynamics_Q
!*
! Compute and print the output of the FEP procedure in Q format
SUBROUTINE AnalyzeFep_Q()
USE Data, ONLY : mappingEnergies
USE Input, ONLY : stateA, coeffs, nFepSteps, mask
USE FreeEnergy, ONLY : ComputeFEPIncrements
IMPLICIT NONE
REAL(8) :: forward(nFepSteps-1), reverse(nFepSteps-1), profile(nFepSteps-1)
REAL(8) :: deltaG_f(nFepSteps), deltaG_r(nFepSteps), deltaG(nFepSteps)
INTEGER :: fepstep
CALL ComputeFEPIncrements(1,nFepSteps,mappingEnergies(:,:,:,1),mask(:,:),forward,reverse,profile)
WRITE(outUnit,'(A43)') "# Part 1: Free energy perturbation summary:"
WRITE(outUnit,*) ""
WRITE(outUnit,'(A)') "# Calculation for full system"
WRITE(outUnit,'(A)') "# lambda(1) dGf sum(dGf) dGr sum(dGr) <dG>"
deltaG_f(:) = 0.0d0
deltaG(:) = 0.0d0
!1st value is 0 in each
DO fepstep = 2, nFepSteps
deltaG_f(fepstep) = forward(fepstep-1)
deltaG(fepstep) = profile(fepstep-1)
ENDDO
!final value is 0
deltaG_r(:) = 0.0d0
DO fepstep = 1, nFepSteps - 1
deltaG_r(fepstep) = reverse(fepstep)
ENDDO
DO fepstep = 1, nFepSteps
WRITE(outUnit,'(2X,F9.6,5F9.3)') coeffs(1,fepstep,stateA), deltaG_f(fepstep), SUM(deltaG_f(1:fepstep)) , &
& deltaG_r(fepstep), SUM(deltaG_r(fepstep:nFepSteps)), &
& SUM(deltaG(1:fepstep))
ENDDO
WRITE(outUnit,*)
ENDSUBROUTINE AnalyzeFep_Q
!*
! Compute and print the output of the FEP/US procedure in Q format
SUBROUTINE AnalyzeFepUs_Q()
USE Data, ONLY : energyGap, mappingEnergies, groundStateEnergy, lambda
USE Input, ONLY : Nbins, minPop, stateA, stateB, stateEnergy, mask
USE FreeEnergy, ONLY : Histogram, ComputeFepProfile, FepUS
IMPLICIT NONE
INTEGER :: bin, step, binPopulations(Nbins,SIZE(energyGap,1)), binIndices(SIZE(energyGap,1),SIZE(energyGap,2))
REAL(8) :: binMidpoints(Nbins)
REAL(8) :: dGg(Nbins,SIZE(energyGap,1)), dGa(Nbins,SIZE(energyGap,1)), dGb(Nbins,SIZE(energyGap,1))
REAL(8) :: binG(Nbins)
REAL(8) :: G_FEP(SIZE(energyGap,1))
INTEGER :: count
WRITE(outUnit,'(A20,F7.2)') "# Min energy-gap is:", MINVAL(energyGap(:,:))
WRITE(outUnit,'(A20,F7.2)') "# Max energy-gap is:", MAXVAL(energyGap(:,:))
WRITE(outUnit,*); WRITE(outUnit,*); WRITE(outUnit,*)
WRITE(outUnit,'(A39)') "# Part 2: Reaction free energy summary:"
WRITE(outUnit,'(A79)') "# Lambda(1) bin Energy gap dGa dGb dGg # pts c1**2 c2**2"
CALL Histogram(energyGap,mask,Nbins,binPopulations,binIndices,binMidpoints)
CALL ComputeFEPProfile(1,SIZE(energyGap,1),mappingEnergies(:,:,:,1),mask(:,:),profile=G_FEP)
! energyGap is the reaction coordinate, nBins is num histogram bins, binPop is histogram values,
! indices is which bin each point is in, binMidpoints is x values
CALL FepUS(mappingEnergies(:,:,:,1),stateEnergy(:,:,stateA,1),G_FEP,binPopulations,binIndices,PMF2Dout=dGa,PMF1D=binG,minPop=minPop)
CALL FepUS(mappingEnergies(:,:,:,1),stateEnergy(:,:,stateB,1),G_FEP,binPopulations,binIndices,PMF2Dout=dGb,PMF1D=binG,minPop=minPop)
CALL FepUS(mappingEnergies(:,:,:,1),groundStateEnergy(:,:), G_FEP,binPopulations,binIndices,PMF2Dout=dGg,PMF1D=binG,minPop=minPop)
DO step = 1, SIZE(energyGap,1)
DO bin = 1, Nbins
IF (binPopulations(bin,step) >= minPop) WRITE(outUnit,'(2X,F9.6,I5,2X,4F9.2,2X,I5)') lambda(step), bin, binMidpoints(bin), &
& dGa(bin,step),dGb(bin,step), dGg(bin,step), binPopulations(bin,step)
ENDDO
ENDDO
WRITE(outUnit,*); WRITE(outUnit,*)
WRITE(outUnit,'(A30)') "# Part 3: Bin-averaged summary"
WRITE(outUnit,'(A63)') "# bin energy gap <dGg> <dGg norm> pts <c1**2> <c2**2> <r_xy>"
DO bin = 1, nBins
count = SUM(binPopulations(bin,:), MASK = binPopulations(bin,:) >= minPop)
IF (SUM(binPopulations(bin,:)) > minPop) WRITE(outUnit,'(I4,1X,3F9.2,2X,I5)') bin, binMidpoints(bin), binG(bin), binG(bin)-MINVAL(binG), count
ENDDO
ENDSUBROUTINE AnalyzeFepUs_Q
END PROGRAM Qfep