forked from cp2k/cp2k
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathfree_energy_methods.F
902 lines (842 loc) · 41.1 KB
/
free_energy_methods.F
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
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
!--------------------------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
! !
! SPDX-License-Identifier: GPL-2.0-or-later !
!--------------------------------------------------------------------------------------------------!
! **************************************************************************************************
!> \brief Methods to perform free energy and free energy derivatives calculations
!> \author Teodoro Laino (01.2007) [tlaino]
! **************************************************************************************************
MODULE free_energy_methods
USE colvar_methods, ONLY: colvar_eval_glob_f
USE cp_log_handling, ONLY: cp_get_default_logger,&
cp_logger_get_default_io_unit,&
cp_logger_type,&
cp_to_string
USE cp_output_handling, ONLY: cp_print_key_finished_output,&
cp_print_key_unit_nr
USE cp_subsys_types, ONLY: cp_subsys_type
USE force_env_types, ONLY: force_env_get,&
force_env_type
USE fparser, ONLY: evalf,&
evalfd,&
finalizef,&
initf,&
parsef
USE free_energy_types, ONLY: free_energy_type,&
ui_var_type
USE input_constants, ONLY: do_fe_ac,&
do_fe_ui
USE input_section_types, ONLY: section_vals_get_subs_vals,&
section_vals_type,&
section_vals_val_get
USE kinds, ONLY: default_path_length,&
default_string_length,&
dp
USE mathlib, ONLY: diamat_all
USE md_environment_types, ONLY: get_md_env,&
md_environment_type
USE memory_utilities, ONLY: reallocate
USE simpar_types, ONLY: simpar_type
USE statistical_methods, ONLY: k_test,&
min_sample_size,&
sw_test,&
vn_test
USE string_utilities, ONLY: compress
#include "../base/base_uses.f90"
IMPLICIT NONE
PRIVATE
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'free_energy_methods'
PUBLIC :: free_energy_evaluate
CONTAINS
! **************************************************************************************************
!> \brief Main driver for free energy calculations
!> In this routine we handle specifically biased MD.
!> \param md_env ...
!> \param converged ...
!> \param fe_section ...
!> \par History
!> Teodoro Laino (01.2007) [tlaino]
! **************************************************************************************************
SUBROUTINE free_energy_evaluate(md_env, converged, fe_section)
TYPE(md_environment_type), POINTER :: md_env
LOGICAL, INTENT(OUT) :: converged
TYPE(section_vals_type), POINTER :: fe_section
CHARACTER(LEN=*), PARAMETER :: routineN = 'free_energy_evaluate'
CHARACTER(LEN=default_path_length) :: coupling_function
CHARACTER(LEN=default_string_length), &
DIMENSION(:), POINTER :: my_par
INTEGER :: handle, ic, icolvar, nforce_eval, &
output_unit, stat_sign_points
INTEGER, POINTER :: istep
REAL(KIND=dp) :: beta, dx, lerr
REAL(KIND=dp), DIMENSION(:), POINTER :: my_val
TYPE(cp_logger_type), POINTER :: logger
TYPE(cp_subsys_type), POINTER :: subsys
TYPE(force_env_type), POINTER :: force_env
TYPE(free_energy_type), POINTER :: fe_env
TYPE(simpar_type), POINTER :: simpar
TYPE(ui_var_type), POINTER :: cv
NULLIFY (force_env, istep, subsys, cv, simpar)
logger => cp_get_default_logger()
CALL timeset(routineN, handle)
converged = .FALSE.
CALL get_md_env(md_env, force_env=force_env, fe_env=fe_env, simpar=simpar, &
itimes=istep)
! Metadynamics is also a free energy calculation but is handled in a different
! module.
IF (.NOT. ASSOCIATED(force_env%meta_env) .AND. ASSOCIATED(fe_env)) THEN
SELECT CASE (fe_env%type)
CASE (do_fe_ui)
! Umbrella Integration..
CALL force_env_get(force_env, subsys=subsys)
fe_env%nr_points = fe_env%nr_points + 1
output_unit = cp_logger_get_default_io_unit(logger)
DO ic = 1, fe_env%ncolvar
cv => fe_env%uivar(ic)
icolvar = cv%icolvar
CALL colvar_eval_glob_f(icolvar, force_env)
CALL reallocate(cv%ss, 1, fe_env%nr_points)
cv%ss(fe_env%nr_points) = subsys%colvar_p(icolvar)%colvar%ss
IF (output_unit > 0) THEN
WRITE (output_unit, *) "COLVAR::", cv%ss(fe_env%nr_points)
END IF
END DO
stat_sign_points = fe_env%nr_points - fe_env%nr_rejected
IF (output_unit > 0) THEN
WRITE (output_unit, *) fe_env%nr_points, stat_sign_points
END IF
! Start statistical analysis when enough CG data points have been collected
IF ((fe_env%conv_par%cg_width*fe_env%conv_par%cg_points <= stat_sign_points) .AND. &
(MOD(stat_sign_points, fe_env%conv_par%cg_width) == 0)) THEN
output_unit = cp_print_key_unit_nr(logger, fe_section, "FREE_ENERGY_INFO", &
extension=".FreeEnergyLog", log_filename=.FALSE.)
CALL print_fe_prolog(output_unit)
! Trend test.. recomputes the number of statistically significant points..
CALL ui_check_trend(fe_env, fe_env%conv_par%test_k, stat_sign_points, output_unit)
stat_sign_points = fe_env%nr_points - fe_env%nr_rejected
! Normality and serial correlation tests..
IF (fe_env%conv_par%cg_width*fe_env%conv_par%cg_points <= stat_sign_points .AND. &
fe_env%conv_par%test_k) THEN
! Statistical tests
CALL ui_check_convergence(fe_env, converged, stat_sign_points, output_unit)
END IF
CALL print_fe_epilog(output_unit)
CALL cp_print_key_finished_output(output_unit, logger, fe_section, "FREE_ENERGY_INFO")
END IF
CASE (do_fe_ac)
CALL initf(2)
! Alchemical Changes
IF (.NOT. ASSOCIATED(force_env%mixed_env)) &
CALL cp_abort(__LOCATION__, &
'ASSERTION (cond) failed at line '//cp_to_string(__LINE__)// &
' Free Energy calculations require the definition of a mixed env!')
my_par => force_env%mixed_env%par
my_val => force_env%mixed_env%val
dx = force_env%mixed_env%dx
lerr = force_env%mixed_env%lerr
coupling_function = force_env%mixed_env%coupling_function
beta = 1/simpar%temp_ext
CALL parsef(1, TRIM(coupling_function), my_par)
nforce_eval = SIZE(force_env%sub_force_env)
CALL dump_ac_info(my_val, my_par, dx, lerr, fe_section, nforce_eval, &
fe_env%covmx, istep, beta)
CALL finalizef()
CASE DEFAULT
! Do Nothing
END SELECT
END IF
CALL timestop(handle)
END SUBROUTINE free_energy_evaluate
! **************************************************************************************************
!> \brief Print prolog of free energy output section
!> \param output_unit which unit to print to
!> \par History
!> Teodoro Laino (02.2007) [tlaino]
! **************************************************************************************************
SUBROUTINE print_fe_prolog(output_unit)
INTEGER, INTENT(IN) :: output_unit
IF (output_unit > 0) THEN
WRITE (output_unit, '(T2,79("*"))')
WRITE (output_unit, '(T30,"FREE ENERGY CALCULATION",/)')
END IF
END SUBROUTINE print_fe_prolog
! **************************************************************************************************
!> \brief Print epilog of free energy output section
!> \param output_unit which unit to print to
!> \par History
!> Teodoro Laino (02.2007) [tlaino]
! **************************************************************************************************
SUBROUTINE print_fe_epilog(output_unit)
INTEGER, INTENT(IN) :: output_unit
IF (output_unit > 0) THEN
WRITE (output_unit, '(T2,79("*"),/)')
END IF
END SUBROUTINE print_fe_epilog
! **************************************************************************************************
!> \brief Test for trend in coarse grained data set
!> \param fe_env ...
!> \param trend_free ...
!> \param nr_points ...
!> \param output_unit which unit to print to
!> \par History
!> Teodoro Laino (01.2007) [tlaino]
! **************************************************************************************************
SUBROUTINE ui_check_trend(fe_env, trend_free, nr_points, output_unit)
TYPE(free_energy_type), POINTER :: fe_env
LOGICAL, INTENT(OUT) :: trend_free
INTEGER, INTENT(IN) :: nr_points, output_unit
CHARACTER(LEN=*), PARAMETER :: routineN = 'ui_check_trend'
INTEGER :: handle, i, ii, j, k, my_reject, ncolvar, &
ng_points, rejected_points
LOGICAL :: test_avg, test_std
REAL(KIND=dp) :: prob, tau, z
REAL(KIND=dp), DIMENSION(:), POINTER :: wrk
CALL timeset(routineN, handle)
trend_free = .FALSE.
test_avg = .TRUE.
test_std = .TRUE.
ncolvar = fe_env%ncolvar
! Number of coarse grained points
IF (output_unit > 0) THEN
WRITE (output_unit, *) nr_points, fe_env%conv_par%cg_width
END IF
ng_points = nr_points/fe_env%conv_par%cg_width
my_reject = 0
! Allocate storage
CALL create_tmp_data(fe_env, wrk, ng_points, ncolvar)
! Compute the Coarse Grained data set using a reverse cumulative strategy
CALL create_csg_data(fe_env, ng_points, output_unit)
! Test on coarse grained average
DO j = 1, ncolvar
ii = 1
DO i = ng_points, 1, -1
wrk(ii) = fe_env%cg_data(i)%avg(j)
ii = ii + 1
END DO
DO i = my_reject + 1, ng_points
IF ((ng_points - my_reject) .LT. min_sample_size) THEN
my_reject = MAX(0, my_reject - 1)
test_avg = .FALSE.
EXIT
END IF
CALL k_test(wrk, my_reject + 1, ng_points, tau, z, prob)
PRINT *, prob, fe_env%conv_par%k_conf_lm
IF (prob < fe_env%conv_par%k_conf_lm) EXIT
my_reject = my_reject + 1
END DO
my_reject = MIN(ng_points, my_reject)
END DO
rejected_points = my_reject*fe_env%conv_par%cg_width
! Print some info
IF (output_unit > 0) THEN
WRITE (output_unit, *) "Kendall trend test (Average)", test_avg, &
"number of points rejected:", rejected_points + fe_env%nr_rejected
WRITE (output_unit, *) "Reject Nr.", my_reject, " coarse grained points testing average"
END IF
! Test on coarse grained covariance matrix
DO j = 1, ncolvar
DO k = j, ncolvar
ii = 1
DO i = ng_points, 1, -1
wrk(ii) = fe_env%cg_data(i)%var(j, k)
ii = ii + 1
END DO
DO i = my_reject + 1, ng_points
IF ((ng_points - my_reject) .LT. min_sample_size) THEN
my_reject = MAX(0, my_reject - 1)
test_std = .FALSE.
EXIT
END IF
CALL k_test(wrk, my_reject + 1, ng_points, tau, z, prob)
PRINT *, prob, fe_env%conv_par%k_conf_lm
IF (prob < fe_env%conv_par%k_conf_lm) EXIT
my_reject = my_reject + 1
END DO
my_reject = MIN(ng_points, my_reject)
END DO
END DO
rejected_points = my_reject*fe_env%conv_par%cg_width
fe_env%nr_rejected = fe_env%nr_rejected + rejected_points
trend_free = test_avg .AND. test_std
! Print some info
IF (output_unit > 0) THEN
WRITE (output_unit, *) "Kendall trend test (Std. Dev.)", test_std, &
"number of points rejected:", fe_env%nr_rejected
WRITE (output_unit, *) "Reject Nr.", my_reject, " coarse grained points testing standard dev."
WRITE (output_unit, *) "Kendall test passed:", trend_free
END IF
! Release storage
CALL destroy_tmp_data(fe_env, wrk, ng_points)
CALL timestop(handle)
END SUBROUTINE ui_check_trend
! **************************************************************************************************
!> \brief Creates temporary data structures
!> \param fe_env ...
!> \param wrk ...
!> \param ng_points ...
!> \param ncolvar ...
!> \par History
!> Teodoro Laino (02.2007) [tlaino]
! **************************************************************************************************
SUBROUTINE create_tmp_data(fe_env, wrk, ng_points, ncolvar)
TYPE(free_energy_type), POINTER :: fe_env
REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: wrk
INTEGER, INTENT(IN) :: ng_points, ncolvar
INTEGER :: i
ALLOCATE (fe_env%cg_data(ng_points))
DO i = 1, ng_points
ALLOCATE (fe_env%cg_data(i)%avg(ncolvar))
ALLOCATE (fe_env%cg_data(i)%var(ncolvar, ncolvar))
END DO
IF (PRESENT(wrk)) THEN
ALLOCATE (wrk(ng_points))
END IF
END SUBROUTINE create_tmp_data
! **************************************************************************************************
!> \brief Destroys temporary data structures
!> \param fe_env ...
!> \param wrk ...
!> \param ng_points ...
!> \par History
!> Teodoro Laino (02.2007) [tlaino]
! **************************************************************************************************
SUBROUTINE destroy_tmp_data(fe_env, wrk, ng_points)
TYPE(free_energy_type), POINTER :: fe_env
REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: wrk
INTEGER, INTENT(IN) :: ng_points
INTEGER :: i
DO i = 1, ng_points
DEALLOCATE (fe_env%cg_data(i)%avg)
DEALLOCATE (fe_env%cg_data(i)%var)
END DO
DEALLOCATE (fe_env%cg_data)
IF (PRESENT(wrk)) THEN
DEALLOCATE (wrk)
END IF
END SUBROUTINE destroy_tmp_data
! **************************************************************************************************
!> \brief Fills in temporary arrays with coarse grained data
!> \param fe_env ...
!> \param ng_points ...
!> \param output_unit which unit to print to
!> \par History
!> Teodoro Laino (02.2007) [tlaino]
! **************************************************************************************************
SUBROUTINE create_csg_data(fe_env, ng_points, output_unit)
TYPE(free_energy_type), POINTER :: fe_env
INTEGER, INTENT(IN) :: ng_points, output_unit
INTEGER :: i, iend, istart
DO i = 1, ng_points
istart = fe_env%nr_points - (i)*fe_env%conv_par%cg_width + 1
iend = fe_env%nr_points - (i - 1)*fe_env%conv_par%cg_width
IF (output_unit > 0) THEN
WRITE (output_unit, *) istart, iend
END IF
CALL eval_cov_matrix(fe_env, cg_index=i, istart=istart, iend=iend, output_unit=output_unit)
END DO
END SUBROUTINE create_csg_data
! **************************************************************************************************
!> \brief Checks Normality of the distribution and Serial Correlation of
!> coarse grained data
!> \param fe_env ...
!> \param test_passed ...
!> \param nr_points ...
!> \param output_unit which unit to print to
!> \par History
!> Teodoro Laino (02.2007) [tlaino]
! **************************************************************************************************
SUBROUTINE ui_check_norm_sc(fe_env, test_passed, nr_points, output_unit)
TYPE(free_energy_type), POINTER :: fe_env
LOGICAL, INTENT(OUT) :: test_passed
INTEGER, INTENT(IN) :: nr_points, output_unit
CHARACTER(LEN=*), PARAMETER :: routineN = 'ui_check_norm_sc'
INTEGER :: handle, ng_points
CALL timeset(routineN, handle)
test_passed = .FALSE.
DO WHILE (fe_env%conv_par%cg_width < fe_env%conv_par%max_cg_width)
ng_points = nr_points/fe_env%conv_par%cg_width
PRINT *, ng_points
IF (ng_points < min_sample_size) EXIT
CALL ui_check_norm_sc_low(fe_env, nr_points, output_unit)
test_passed = fe_env%conv_par%test_vn .AND. fe_env%conv_par%test_sw
IF (test_passed) EXIT
fe_env%conv_par%cg_width = fe_env%conv_par%cg_width + 1
IF (output_unit > 0) THEN
WRITE (output_unit, *) "New coarse grained width:", fe_env%conv_par%cg_width
END IF
END DO
IF (fe_env%conv_par%cg_width == fe_env%conv_par%max_cg_width .AND. (.NOT. (test_passed))) THEN
CALL ui_check_norm_sc_low(fe_env, nr_points, output_unit)
test_passed = fe_env%conv_par%test_vn .AND. fe_env%conv_par%test_sw
END IF
CALL timestop(handle)
END SUBROUTINE ui_check_norm_sc
! **************************************************************************************************
!> \brief Checks Normality of the distribution and Serial Correlation of
!> coarse grained data - Low Level routine
!> \param fe_env ...
!> \param nr_points ...
!> \param output_unit which unit to print to
!> \par History
!> Teodoro Laino (02.2007) [tlaino]
! **************************************************************************************************
SUBROUTINE ui_check_norm_sc_low(fe_env, nr_points, output_unit)
TYPE(free_energy_type), POINTER :: fe_env
INTEGER, INTENT(IN) :: nr_points, output_unit
CHARACTER(LEN=*), PARAMETER :: routineN = 'ui_check_norm_sc_low'
INTEGER :: handle, i, j, k, ncolvar, ng_points
LOGICAL :: avg_test_passed, sdv_test_passed
REAL(KIND=dp) :: prob, pw, r, u, w
REAL(KIND=dp), DIMENSION(:), POINTER :: wrk
CALL timeset(routineN, handle)
ncolvar = fe_env%ncolvar
! Compute the Coarse Grained data set using a reverse cumulative strategy
fe_env%conv_par%test_sw = .FALSE.
fe_env%conv_par%test_vn = .FALSE.
! Number of coarse grained points
avg_test_passed = .TRUE.
sdv_test_passed = .TRUE.
ng_points = nr_points/fe_env%conv_par%cg_width
CALL create_tmp_data(fe_env, wrk, ng_points, ncolvar)
CALL create_csg_data(fe_env, ng_points, output_unit)
! Testing Averages
DO j = 1, ncolvar
DO i = 1, ng_points
wrk(i) = fe_env%cg_data(i)%avg(j)
END DO
! Test of Shapiro - Wilks for normality
! - Average
CALL sw_test(wrk, ng_points, w, pw)
PRINT *, 1.0_dp - pw, fe_env%conv_par%sw_conf_lm
avg_test_passed = (1.0_dp - pw) <= fe_env%conv_par%sw_conf_lm
fe_env%conv_par%test_sw = avg_test_passed
IF (output_unit > 0) THEN
WRITE (output_unit, *) "Shapiro-Wilks normality test (Avg)", avg_test_passed
END IF
! Test of von Neumann for serial correlation
! - Average
CALL vn_test(wrk, ng_points, r, u, prob)
PRINT *, prob, fe_env%conv_par%vn_conf_lm
avg_test_passed = prob <= fe_env%conv_par%vn_conf_lm
fe_env%conv_par%test_vn = avg_test_passed
IF (output_unit > 0) THEN
WRITE (output_unit, *) "von Neumann serial correlation test (Avg)", avg_test_passed
END IF
END DO
! If tests on average are ok let's proceed with Standard Deviation
IF (fe_env%conv_par%test_vn .AND. fe_env%conv_par%test_sw) THEN
! Testing Standard Deviations
DO j = 1, ncolvar
DO k = j, ncolvar
DO i = 1, ng_points
wrk(i) = fe_env%cg_data(i)%var(j, k)
END DO
! Test of Shapiro - Wilks for normality
! - Standard Deviation
CALL sw_test(wrk, ng_points, w, pw)
PRINT *, 1.0_dp - pw, fe_env%conv_par%sw_conf_lm
sdv_test_passed = (1.0_dp - pw) <= fe_env%conv_par%sw_conf_lm
fe_env%conv_par%test_sw = fe_env%conv_par%test_sw .AND. sdv_test_passed
IF (output_unit > 0) THEN
WRITE (output_unit, *) "Shapiro-Wilks normality test (Std. Dev.)", sdv_test_passed
END IF
! Test of von Neumann for serial correlation
! - Standard Deviation
CALL vn_test(wrk, ng_points, r, u, prob)
PRINT *, prob, fe_env%conv_par%vn_conf_lm
sdv_test_passed = prob <= fe_env%conv_par%vn_conf_lm
fe_env%conv_par%test_vn = fe_env%conv_par%test_vn .AND. sdv_test_passed
IF (output_unit > 0) THEN
WRITE (output_unit, *) "von Neumann serial correlation test (Std. Dev.)", sdv_test_passed
END IF
END DO
END DO
CALL destroy_tmp_data(fe_env, wrk, ng_points)
ELSE
CALL destroy_tmp_data(fe_env, wrk, ng_points)
END IF
CALL timestop(handle)
END SUBROUTINE ui_check_norm_sc_low
! **************************************************************************************************
!> \brief Convergence criteria (Error on average and covariance matrix)
!> for free energy method
!> \param fe_env ...
!> \param converged ...
!> \param nr_points ...
!> \param output_unit which unit to print to
!> \par History
!> Teodoro Laino (01.2007) [tlaino]
! **************************************************************************************************
SUBROUTINE ui_check_convergence(fe_env, converged, nr_points, output_unit)
TYPE(free_energy_type), POINTER :: fe_env
LOGICAL, INTENT(OUT) :: converged
INTEGER, INTENT(IN) :: nr_points, output_unit
CHARACTER(LEN=*), PARAMETER :: routineN = 'ui_check_convergence'
INTEGER :: handle, i, ic, ncolvar, ng_points
LOGICAL :: test_passed
REAL(KIND=dp) :: max_error_avg, max_error_std
REAL(KIND=dp), DIMENSION(:), POINTER :: avg_std, avgmx
REAL(KIND=dp), DIMENSION(:, :), POINTER :: cov_std, covmx
CALL timeset(routineN, handle)
converged = .FALSE.
ncolvar = fe_env%ncolvar
NULLIFY (avgmx, avg_std, covmx, cov_std)
CALL ui_check_norm_sc(fe_env, test_passed, nr_points, output_unit)
IF (test_passed) THEN
ng_points = nr_points/fe_env%conv_par%cg_width
! We can finally compute the error on average and covariance matrix
! and check if we converged..
CALL create_tmp_data(fe_env, ng_points=ng_points, ncolvar=ncolvar)
CALL create_csg_data(fe_env, ng_points, output_unit)
ALLOCATE (covmx(ncolvar, ncolvar))
ALLOCATE (avgmx(ncolvar))
ALLOCATE (cov_std(ncolvar*(ncolvar + 1)/2, ncolvar*(ncolvar + 1)/2))
ALLOCATE (avg_std(ncolvar))
covmx = 0.0_dp
avgmx = 0.0_dp
DO i = 1, ng_points
covmx = covmx + fe_env%cg_data(i)%var
avgmx = avgmx + fe_env%cg_data(i)%avg
END DO
covmx = covmx/REAL(ng_points, KIND=dp)
avgmx = avgmx/REAL(ng_points, KIND=dp)
! Compute errors on average and standard deviation
CALL compute_avg_std_errors(fe_env, ncolvar, avgmx, covmx, avg_std, cov_std)
IF (output_unit > 0) THEN
WRITE (output_unit, *) "pippo", avgmx, covmx
WRITE (output_unit, *) "pippo", avg_std, cov_std
END IF
! Convergence of the averages
max_error_avg = SQRT(MAXVAL(ABS(avg_std))/REAL(ng_points, KIND=dp))/MINVAL(avgmx)
max_error_std = SQRT(MAXVAL(ABS(cov_std))/REAL(ng_points, KIND=dp))/MINVAL(covmx)
IF (max_error_avg <= fe_env%conv_par%eps_conv .AND. &
max_error_std <= fe_env%conv_par%eps_conv) converged = .TRUE.
IF (output_unit > 0) THEN
WRITE (output_unit, '(/,T2,"CG SAMPLING LENGTH = ",I7,20X,"REQUESTED ACCURACY = ",E12.6)') ng_points, &
fe_env%conv_par%eps_conv
WRITE (output_unit, '(T50,"PRESENT ACCURACY AVG= ",E12.6)') max_error_avg
WRITE (output_unit, '(T50,"PRESENT ACCURACY STD= ",E12.6)') max_error_std
WRITE (output_unit, '(T50,"CONVERGED FE-DER = ",L12)') converged
WRITE (output_unit, '(/,T33, "COVARIANCE MATRIX")')
WRITE (output_unit, '(T8,'//cp_to_string(ncolvar)//'(3X,I7,6X))') (ic, ic=1, ncolvar)
DO ic = 1, ncolvar
WRITE (output_unit, '(T2,I6,'//cp_to_string(ncolvar)//'(3X,E12.6,1X))') ic, covmx(ic, :)
END DO
WRITE (output_unit, '(T33, "ERROR OF COVARIANCE MATRIX")')
WRITE (output_unit, '(T8,'//cp_to_string(ncolvar)//'(3X,I7,6X))') (ic, ic=1, ncolvar)
DO ic = 1, ncolvar
WRITE (output_unit, '(T2,I6,'//cp_to_string(ncolvar)//'(3X,E12.6,1X))') ic, cov_std(ic, :)
END DO
WRITE (output_unit, '(/,T2,"COLVAR Nr.",18X,13X,"AVERAGE",13X,"STANDARD DEVIATION")')
WRITE (output_unit, '(T2,"CV",I8,21X,7X,E12.6,14X,E12.6)') &
(ic, avgmx(ic), SQRT(ABS(avg_std(ic))), ic=1, ncolvar)
END IF
CALL destroy_tmp_data(fe_env, ng_points=ng_points)
DEALLOCATE (covmx)
DEALLOCATE (avgmx)
DEALLOCATE (cov_std)
DEALLOCATE (avg_std)
END IF
CALL timestop(handle)
END SUBROUTINE ui_check_convergence
! **************************************************************************************************
!> \brief Computes the errors on averages and standard deviations for a
!> correlation-independent coarse grained data set
!> \param fe_env ...
!> \param ncolvar ...
!> \param avgmx ...
!> \param covmx ...
!> \param avg_std ...
!> \param cov_std ...
!> \par History
!> Teodoro Laino (02.2007) [tlaino]
! **************************************************************************************************
SUBROUTINE compute_avg_std_errors(fe_env, ncolvar, avgmx, covmx, avg_std, cov_std)
TYPE(free_energy_type), POINTER :: fe_env
INTEGER, INTENT(IN) :: ncolvar
REAL(KIND=dp), DIMENSION(:), POINTER :: avgmx
REAL(KIND=dp), DIMENSION(:, :), POINTER :: covmx
REAL(KIND=dp), DIMENSION(:), POINTER :: avg_std
REAL(KIND=dp), DIMENSION(:, :), POINTER :: cov_std
INTEGER :: i, ind, j, k, nvar
REAL(KIND=dp) :: fac
REAL(KIND=dp), DIMENSION(:), POINTER :: awrk, eig, tmp
REAL(KIND=dp), DIMENSION(:, :), POINTER :: wrk
! Averages
nvar = ncolvar
ALLOCATE (wrk(nvar, nvar))
ALLOCATE (eig(nvar))
fac = REAL(SIZE(fe_env%cg_data), KIND=dp)
wrk = 0.0_dp
eig = 0.0_dp
DO k = 1, SIZE(fe_env%cg_data)
DO j = 1, nvar
DO i = j, nvar
wrk(i, j) = wrk(i, j) + fe_env%cg_data(k)%avg(i)*fe_env%cg_data(k)%avg(j)
END DO
END DO
END DO
DO j = 1, nvar
DO i = j, nvar
wrk(i, j) = wrk(i, j) - avgmx(i)*avgmx(j)*fac
wrk(j, i) = wrk(i, j)
END DO
END DO
wrk = wrk/(fac - 1.0_dp)
! Diagonalize the covariance matrix and check for the maximum error
CALL diamat_all(wrk, eig)
DO i = 1, nvar
avg_std(i) = eig(i)
END DO
DEALLOCATE (wrk)
DEALLOCATE (eig)
! Standard Deviations
nvar = ncolvar*(ncolvar + 1)/2
ALLOCATE (wrk(nvar, nvar))
ALLOCATE (eig(nvar))
ALLOCATE (awrk(nvar))
ALLOCATE (tmp(nvar))
wrk = 0.0_dp
eig = 0.0_dp
ind = 0
DO i = 1, ncolvar
DO j = i, ncolvar
ind = ind + 1
awrk(ind) = covmx(i, j)
END DO
END DO
DO k = 1, SIZE(fe_env%cg_data)
ind = 0
DO i = 1, ncolvar
DO j = i, ncolvar
ind = ind + 1
tmp(ind) = fe_env%cg_data(k)%var(i, j)
END DO
END DO
DO i = 1, nvar
DO j = i, nvar
wrk(i, j) = wrk(i, j) + tmp(i)*tmp(j) - awrk(i)*awrk(j)
END DO
END DO
END DO
DO i = 1, nvar
DO j = i, nvar
wrk(i, j) = wrk(i, j) - fac*awrk(i)*awrk(j)
wrk(j, i) = wrk(i, j)
END DO
END DO
wrk = wrk/(fac - 1.0_dp)
! Diagonalize the covariance matrix and check for the maximum error
CALL diamat_all(wrk, eig)
ind = 0
DO i = 1, ncolvar
DO j = i, ncolvar
ind = ind + 1
cov_std(i, j) = eig(ind)
cov_std(j, i) = cov_std(i, j)
END DO
END DO
DEALLOCATE (wrk)
DEALLOCATE (eig)
DEALLOCATE (awrk)
DEALLOCATE (tmp)
END SUBROUTINE compute_avg_std_errors
! **************************************************************************************************
!> \brief Computes the covariance matrix
!> \param fe_env ...
!> \param cg_index ...
!> \param istart ...
!> \param iend ...
!> \param output_unit which unit to print to
!> \param covmx ...
!> \param avgs ...
!> \par History
!> Teodoro Laino (01.2007) [tlaino]
! **************************************************************************************************
SUBROUTINE eval_cov_matrix(fe_env, cg_index, istart, iend, output_unit, covmx, avgs)
TYPE(free_energy_type), POINTER :: fe_env
INTEGER, INTENT(IN) :: cg_index, istart, iend, output_unit
REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER :: covmx
REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: avgs
CHARACTER(LEN=*), PARAMETER :: routineN = 'eval_cov_matrix'
INTEGER :: handle, ic, jc, jstep, ncolvar, nlength
REAL(KIND=dp) :: tmp_ic, tmp_jc
TYPE(ui_var_type), POINTER :: cv
CALL timeset(routineN, handle)
ncolvar = fe_env%ncolvar
nlength = iend - istart + 1
fe_env%cg_data(cg_index)%avg = 0.0_dp
fe_env%cg_data(cg_index)%var = 0.0_dp
IF (nlength > 1) THEN
! Update the info on averages and variances
DO jstep = istart, iend
DO ic = 1, ncolvar
cv => fe_env%uivar(ic)
tmp_ic = cv%ss(jstep)
fe_env%cg_data(cg_index)%avg(ic) = fe_env%cg_data(cg_index)%avg(ic) + tmp_ic
END DO
DO ic = 1, ncolvar
cv => fe_env%uivar(ic)
tmp_ic = cv%ss(jstep)
DO jc = 1, ic
cv => fe_env%uivar(jc)
tmp_jc = cv%ss(jstep)
fe_env%cg_data(cg_index)%var(jc, ic) = fe_env%cg_data(cg_index)%var(jc, ic) + tmp_ic*tmp_jc
END DO
END DO
END DO
! Normalized the variances and the averages
! Unbiased estimator
fe_env%cg_data(cg_index)%var = fe_env%cg_data(cg_index)%var/REAL(nlength - 1, KIND=dp)
fe_env%cg_data(cg_index)%avg = fe_env%cg_data(cg_index)%avg/REAL(nlength, KIND=dp)
! Compute the covariance matrix
DO ic = 1, ncolvar
tmp_ic = fe_env%cg_data(cg_index)%avg(ic)
DO jc = 1, ic
tmp_jc = fe_env%cg_data(cg_index)%avg(jc)*REAL(nlength, KIND=dp)/REAL(nlength - 1, KIND=dp)
fe_env%cg_data(cg_index)%var(jc, ic) = fe_env%cg_data(cg_index)%var(jc, ic) - tmp_ic*tmp_jc
fe_env%cg_data(cg_index)%var(ic, jc) = fe_env%cg_data(cg_index)%var(jc, ic)
END DO
END DO
IF (output_unit > 0) THEN
WRITE (output_unit, *) "eval_cov_matrix", istart, iend, fe_env%cg_data(cg_index)%avg, fe_env%cg_data(cg_index)%var
END IF
IF (PRESENT(covmx)) covmx = fe_env%cg_data(cg_index)%var
IF (PRESENT(avgs)) avgs = fe_env%cg_data(cg_index)%avg
END IF
CALL timestop(handle)
END SUBROUTINE eval_cov_matrix
! **************************************************************************************************
!> \brief Dumps information when performing an alchemical change run
!> \param my_val ...
!> \param my_par ...
!> \param dx ...
!> \param lerr ...
!> \param fe_section ...
!> \param nforce_eval ...
!> \param cum_res ...
!> \param istep ...
!> \param beta ...
!> \author Teodoro Laino - University of Zurich [tlaino] - 05.2007
! **************************************************************************************************
SUBROUTINE dump_ac_info(my_val, my_par, dx, lerr, fe_section, nforce_eval, cum_res, &
istep, beta)
REAL(KIND=dp), DIMENSION(:), POINTER :: my_val
CHARACTER(LEN=default_string_length), &
DIMENSION(:), POINTER :: my_par
REAL(KIND=dp), INTENT(IN) :: dx, lerr
TYPE(section_vals_type), POINTER :: fe_section
INTEGER, INTENT(IN) :: nforce_eval
REAL(KIND=dp), DIMENSION(:, :), POINTER :: cum_res
INTEGER, POINTER :: istep
REAL(KIND=dp), INTENT(IN) :: beta
CHARACTER(LEN=default_path_length) :: coupling_function
CHARACTER(LEN=default_string_length) :: def_error, par, this_error
INTEGER :: i, iforce_eval, ipar, isize, iw, j, &
NEquilStep
REAL(KIND=dp) :: avg_BP, avg_DET, avg_DUE, d_ene_w, dedf, &
ene_w, err, Err_DET, Err_DUE, std_DET, &
std_DUE, tmp, tmp2, wfac
TYPE(cp_logger_type), POINTER :: logger
TYPE(section_vals_type), POINTER :: alch_section
logger => cp_get_default_logger()
alch_section => section_vals_get_subs_vals(fe_section, "ALCHEMICAL_CHANGE")
CALL section_vals_val_get(alch_section, "PARAMETER", c_val=par)
DO i = 1, SIZE(my_par)
IF (my_par(i) == par) EXIT
END DO
CPASSERT(i <= SIZE(my_par))
ipar = i
dedf = evalfd(1, ipar, my_val, dx, err)
IF (ABS(err) > lerr) THEN
WRITE (this_error, "(A,G12.6,A)") "(", err, ")"
WRITE (def_error, "(A,G12.6,A)") "(", lerr, ")"
CALL compress(this_error, .TRUE.)
CALL compress(def_error, .TRUE.)
CALL cp_warn(__LOCATION__, &
'ASSERTION (cond) failed at line '//cp_to_string(__LINE__)// &
' Error '//TRIM(this_error)//' in computing numerical derivatives larger then'// &
TRIM(def_error)//' .')
END IF
! We must print now the energy of the biased system, the weigthing energy
! and the derivative w.r.t.the coupling parameter of the biased energy
! Retrieve the expression of the weighting function:
CALL section_vals_val_get(alch_section, "WEIGHTING_FUNCTION", c_val=coupling_function)
CALL compress(coupling_function, full=.TRUE.)
CALL parsef(2, TRIM(coupling_function), my_par)
ene_w = evalf(2, my_val)
d_ene_w = evalfd(2, ipar, my_val, dx, err)
IF (ABS(err) > lerr) THEN
WRITE (this_error, "(A,G12.6,A)") "(", err, ")"
WRITE (def_error, "(A,G12.6,A)") "(", lerr, ")"
CALL compress(this_error, .TRUE.)
CALL compress(def_error, .TRUE.)
CALL cp_warn(__LOCATION__, &
'ASSERTION (cond) failed at line '//cp_to_string(__LINE__)// &
' Error '//TRIM(this_error)//' in computing numerical derivatives larger then'// &
TRIM(def_error)//' .')
END IF
CALL section_vals_val_get(alch_section, "NEQUIL_STEPS", i_val=NEquilStep)
! Store results
IF (istep > NEquilStep) THEN
isize = SIZE(cum_res, 2) + 1
CALL reallocate(cum_res, 1, 3, 1, isize)
cum_res(1, isize) = dedf
cum_res(2, isize) = dedf - d_ene_w
cum_res(3, isize) = ene_w
! Compute derivative of biased and total energy
! Total Free Energy
avg_DET = SUM(cum_res(1, 1:isize))/REAL(isize, KIND=dp)
std_DET = SUM(cum_res(1, 1:isize)**2)/REAL(isize, KIND=dp)
! Unbiased Free Energy
avg_BP = SUM(cum_res(3, 1:isize))/REAL(isize, KIND=dp)
wfac = 0.0_dp
DO j = 1, isize
wfac = wfac + EXP(beta*(cum_res(3, j) - avg_BP))
END DO
avg_DUE = 0.0_dp
std_DUE = 0.0_dp
DO j = 1, isize
tmp = cum_res(2, j)
tmp2 = EXP(beta*(cum_res(3, j) - avg_BP))/wfac
avg_DUE = avg_DUE + tmp*tmp2
std_DUE = std_DUE + tmp**2*tmp2
END DO
IF (isize > 1) THEN
Err_DUE = SQRT(std_DUE - avg_DUE**2)/SQRT(REAL(isize - 1, KIND=dp))
Err_DET = SQRT(std_DET - avg_DET**2)/SQRT(REAL(isize - 1, KIND=dp))
END IF
! Print info
iw = cp_print_key_unit_nr(logger, fe_section, "FREE_ENERGY_INFO", &
extension=".free_energy")
IF (iw > 0) THEN
WRITE (iw, '(T2,79("-"),T37," oOo ")')
DO iforce_eval = 1, nforce_eval
WRITE (iw, '(T2,"ALCHEMICAL CHANGE| FORCE_EVAL Nr.",I5,T48,"ENERGY (Hartree)= ",F15.9)') &
iforce_eval, my_val(iforce_eval)
END DO
WRITE (iw, '(T2,"ALCHEMICAL CHANGE| DERIVATIVE OF TOTAL ENERGY [ PARAMETER (",A,") ]",T66,F15.9)') &
TRIM(par), dedf
WRITE (iw, '(T2,"ALCHEMICAL CHANGE| DERIVATIVE OF BIASED ENERGY [ PARAMETER (",A,") ]",T66,F15.9)') &
TRIM(par), dedf - d_ene_w
WRITE (iw, '(T2,"ALCHEMICAL CHANGE| BIASING UMBRELLA POTENTIAL ",T66,F15.9)') &
ene_w
IF (isize > 1) THEN
WRITE (iw, '(/,T2,"ALCHEMICAL CHANGE| DERIVATIVE TOTAL FREE ENERGY ",T50,F15.9,1X,"+/-",1X,F11.9)') &
avg_DET, Err_DET
WRITE (iw, '(T2,"ALCHEMICAL CHANGE| DERIVATIVE UNBIASED FREE ENERGY ",T50,F15.9,1X,"+/-",1X,F11.9)') &
avg_DUE, Err_DUE
ELSE
WRITE (iw, '(/,T2,"ALCHEMICAL CHANGE| DERIVATIVE TOTAL FREE ENERGY ",T50,F15.9,1X,"+/-",1X,T76,A)') &
avg_DET, "UNDEF"
WRITE (iw, '(T2,"ALCHEMICAL CHANGE| DERIVATIVE UNBIASED FREE ENERGY ",T50,F15.9,1X,"+/-",1X,T76,A)') &
avg_DUE, "UNDEF"
END IF
WRITE (iw, '(T2,79("-"))')
END IF
END IF
CALL cp_print_key_finished_output(iw, logger, fe_section, "FREE_ENERGY_INFO")
END SUBROUTINE dump_ac_info
END MODULE free_energy_methods