-
Notifications
You must be signed in to change notification settings - Fork 1
/
metadynamics.F
1069 lines (949 loc) · 47.3 KB
/
metadynamics.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
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
!--------------------------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
! !
! SPDX-License-Identifier: GPL-2.0-or-later !
!--------------------------------------------------------------------------------------------------!
! **************************************************************************************************
!> \brief Performs the metadynamics calculation
!> \par History
!> 01.2005 created [fawzi and ale]
!> 11.2007 Teodoro Laino [tlaino] - University of Zurich
! **************************************************************************************************
MODULE metadynamics
USE cp_log_handling, ONLY: cp_logger_type, cp_get_default_logger
USE bibliography, ONLY: VandenCic2006
#if defined (__PLUMED2)
USE ISO_C_BINDING, ONLY: C_INT, C_DOUBLE, C_CHAR
USE cell_types, ONLY: cell_type, &
pbc_cp2k_plumed_getset_cell
#else
USE cell_types, ONLY: cell_type
#endif
USE colvar_methods, ONLY: colvar_eval_glob_f
USE colvar_types, ONLY: colvar_p_type, &
torsion_colvar_id
USE constraint_fxd, ONLY: fix_atom_control
USE cp_output_handling, ONLY: cp_add_iter_level, &
cp_iterate, &
cp_print_key_finished_output, &
cp_print_key_unit_nr, &
cp_rm_iter_level
USE cp_subsys_types, ONLY: cp_subsys_get, &
cp_subsys_type
USE force_env_types, ONLY: force_env_get, &
force_env_type
USE input_constants, ONLY: do_wall_m, &
do_wall_p, &
do_wall_reflective
USE input_section_types, ONLY: section_vals_get, &
section_vals_get_subs_vals, &
section_vals_type, &
section_vals_val_get
USE kinds, ONLY: dp
USE message_passing, ONLY: cp2k_is_parallel
USE metadynamics_types, ONLY: hills_env_type, &
meta_env_type, &
metavar_type, &
multiple_walkers_type
USE metadynamics_utils, ONLY: add_hill_single, &
get_meta_iter_level, &
meta_walls, &
restart_hills, &
synchronize_multiple_walkers
USE particle_list_types, ONLY: particle_list_type
#if defined (__PLUMED2)
USE physcon, ONLY: angstrom, &
boltzmann, &
femtoseconds, &
joule, &
kelvin, kjmol, &
picoseconds
#else
USE physcon, ONLY: boltzmann, &
femtoseconds, &
joule, &
kelvin
#endif
USE reference_manager, ONLY: cite_reference
USE simpar_types, ONLY: simpar_type
#include "./base/base_uses.f90"
IMPLICIT NONE
PRIVATE
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'metadynamics'
PUBLIC :: metadyn_forces, metadyn_integrator
PUBLIC :: metadyn_velocities_colvar, metadyn_write_colvar
PUBLIC :: metadyn_initialise_plumed, metadyn_finalise_plumed
#if defined (__PLUMED2)
INTERFACE
FUNCTION plumed_installed() RESULT(res) BIND(C, name="plumed_installed")
IMPORT :: C_INT
INTEGER(KIND=C_INT) :: res
END FUNCTION plumed_installed
SUBROUTINE plumed_gcreate() BIND(C, name="plumed_gcreate")
END SUBROUTINE plumed_gcreate
SUBROUTINE plumed_gfinalize() BIND(C, name="plumed_gfinalize")
END SUBROUTINE plumed_gfinalize
SUBROUTINE plumed_gcmd_int(key, value) BIND(C, name="plumed_gcmd")
IMPORT :: C_CHAR, C_INT
CHARACTER(KIND=C_CHAR), DIMENSION(*) :: key
INTEGER(KIND=C_INT) :: value
END SUBROUTINE plumed_gcmd_int
SUBROUTINE plumed_gcmd_double(key, value) BIND(C, name="plumed_gcmd")
IMPORT :: C_CHAR, C_DOUBLE
CHARACTER(KIND=C_CHAR), DIMENSION(*) :: key
REAL(KIND=C_DOUBLE) :: value
END SUBROUTINE plumed_gcmd_double
SUBROUTINE plumed_gcmd_doubles(key, value) BIND(C, name="plumed_gcmd")
IMPORT :: C_CHAR, C_DOUBLE
CHARACTER(KIND=C_CHAR), DIMENSION(*) :: key
REAL(KIND=C_DOUBLE), DIMENSION(*) :: value
END SUBROUTINE plumed_gcmd_doubles
SUBROUTINE plumed_gcmd_char(key, value) BIND(C, name="plumed_gcmd")
IMPORT :: C_CHAR
CHARACTER(KIND=C_CHAR), DIMENSION(*) :: key, value
END SUBROUTINE plumed_gcmd_char
END INTERFACE
#endif
CONTAINS
! **************************************************************************************************
!> \brief ...
!> \param force_env ...
!> \param simpar ...
!> \param itimes ...
! **************************************************************************************************
SUBROUTINE metadyn_initialise_plumed(force_env, simpar, itimes)
TYPE(force_env_type), POINTER :: force_env
TYPE(simpar_type), POINTER :: simpar
INTEGER, INTENT(IN) :: itimes
CHARACTER(LEN=*), PARAMETER :: routineN = 'metadyn_initialise_plumed'
INTEGER :: handle
#if defined (__PLUMED2)
INTEGER :: natom_plumed
REAL(KIND=dp) :: timestep_plumed
TYPE(cell_type), POINTER :: cell
TYPE(cp_subsys_type), POINTER :: subsys
#endif
#if defined (__PLUMED2)
INTEGER :: plumedavailable
#endif
CALL timeset(routineN, handle)
CPASSERT(ASSOCIATED(force_env))
CPASSERT(ASSOCIATED(simpar))
#if defined (__PLUMED2)
NULLIFY (cell, subsys)
CALL force_env_get(force_env, subsys=subsys, cell=cell)
CALL pbc_cp2k_plumed_getset_cell(cell, set=.TRUE.) !Store the cell pointer for later use.
timestep_plumed = simpar%dt
natom_plumed = subsys%particles%n_els
#else
MARK_USED(force_env)
MARK_USED(simpar)
#endif
MARK_USED(itimes)
#if defined (__PLUMED2)
plumedavailable = plumed_installed()
IF (plumedavailable <= 0) THEN
CPABORT("NO PLUMED IN YOUR LD_LIBRARY_PATH?")
ELSE
CALL plumed_gcreate()
IF (cp2k_is_parallel) CALL plumed_gcmd_int("setMPIFComm"//CHAR(0), force_env%para_env%get_handle())
CALL plumed_gcmd_int("setRealPrecision"//CHAR(0), 8)
CALL plumed_gcmd_double("setMDEnergyUnits"//CHAR(0), kjmol) ! Hartree to kjoule/mol
CALL plumed_gcmd_double("setMDLengthUnits"//CHAR(0), angstrom*0.1_dp) ! Bohr to nm
CALL plumed_gcmd_double("setMDTimeUnits"//CHAR(0), picoseconds) ! atomic units to ps
CALL plumed_gcmd_char("setPlumedDat"//CHAR(0), force_env%meta_env%plumed_input_file//CHAR(0))
CALL plumed_gcmd_char("setLogFile"//CHAR(0), "PLUMED.OUT"//CHAR(0))
CALL plumed_gcmd_int("setNatoms"//CHAR(0), natom_plumed)
CALL plumed_gcmd_char("setMDEngine"//CHAR(0), "cp2k"//CHAR(0))
CALL plumed_gcmd_double("setTimestep"//CHAR(0), timestep_plumed)
CALL plumed_gcmd_int("init"//CHAR(0), 0)
END IF
#endif
#if !defined (__PLUMED2)
CALL cp_abort(__LOCATION__, &
"Requested to use plumed for metadynamics, but cp2k was"// &
" not compiled with plumed support.")
#endif
CALL timestop(handle)
END SUBROUTINE
! **************************************************************************************************
!> \brief makes sure PLUMED is shut down cleanly
! **************************************************************************************************
SUBROUTINE metadyn_finalise_plumed()
CHARACTER(LEN=*), PARAMETER :: routineN = 'metadyn_finalise_plumed'
INTEGER :: handle
CALL timeset(routineN, handle)
#if defined (__PLUMED2)
CALL plumed_gfinalize()
#endif
CALL timestop(handle)
END SUBROUTINE
! **************************************************************************************************
!> \brief General driver for applying metadynamics
!> \param force_env ...
!> \param itimes ...
!> \param vel ...
!> \param rand ...
!> \date 01.2009
!> \par History
!> 01.2009 created
!> \author Teodoro Laino
! **************************************************************************************************
SUBROUTINE metadyn_integrator(force_env, itimes, vel, rand)
TYPE(force_env_type), POINTER :: force_env
INTEGER, INTENT(IN) :: itimes
REAL(KIND=dp), DIMENSION(:, :), &
INTENT(INOUT), OPTIONAL :: vel
REAL(KIND=dp), DIMENSION(:), OPTIONAL, &
POINTER :: rand
CHARACTER(len=*), PARAMETER :: routineN = 'metadyn_integrator'
INTEGER :: handle, plumed_itimes
#if defined (__PLUMED2)
INTEGER :: i_part, natom_plumed
TYPE(cell_type), POINTER :: cell
TYPE(cp_subsys_type), POINTER :: subsys
#endif
#if defined (__PLUMED2)
TYPE(meta_env_type), POINTER :: meta_env
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: mass_plumed
REAL(KIND=dp), DIMENSION(3, 3) :: plu_vir, celltbt
REAL(KIND=dp) :: stpcfg
REAL(KIND=dp), ALLOCATABLE, &
DIMENSION(:) :: pos_plumed_x, pos_plumed_y, pos_plumed_z
REAL(KIND=dp), ALLOCATABLE, &
DIMENSION(:) :: force_plumed_x, force_plumed_y, force_plumed_z
#endif
CALL timeset(routineN, handle)
! Apply Metadynamics
IF (ASSOCIATED(force_env%meta_env)) THEN
IF (force_env%meta_env%use_plumed .EQV. .TRUE.) THEN
plumed_itimes = itimes
#if defined (__PLUMED2)
CALL force_env_get(force_env, meta_env=meta_env, subsys=subsys, cell=cell)
natom_plumed = subsys%particles%n_els
ALLOCATE (pos_plumed_x(natom_plumed))
ALLOCATE (pos_plumed_y(natom_plumed))
ALLOCATE (pos_plumed_z(natom_plumed))
ALLOCATE (force_plumed_x(natom_plumed))
ALLOCATE (force_plumed_y(natom_plumed))
ALLOCATE (force_plumed_z(natom_plumed))
ALLOCATE (mass_plumed(natom_plumed))
force_plumed_x(:) = 0.0d0
force_plumed_y(:) = 0.0d0
force_plumed_z(:) = 0.0d0
DO i_part = 1, natom_plumed
pos_plumed_x(i_part) = subsys%particles%els(i_part)%r(1)
pos_plumed_y(i_part) = subsys%particles%els(i_part)%r(2)
pos_plumed_z(i_part) = subsys%particles%els(i_part)%r(3)
mass_plumed(i_part) = subsys%particles%els(i_part)%atomic_kind%mass
END DO
! stupid cell conversion, to be fixed
celltbt(1, 1) = cell%hmat(1, 1)
celltbt(1, 2) = cell%hmat(1, 2)
celltbt(1, 3) = cell%hmat(1, 3)
celltbt(2, 1) = cell%hmat(2, 1)
celltbt(2, 2) = cell%hmat(2, 2)
celltbt(2, 3) = cell%hmat(2, 3)
celltbt(3, 1) = cell%hmat(3, 1)
celltbt(3, 2) = cell%hmat(3, 2)
celltbt(3, 3) = cell%hmat(3, 3)
! virial
plu_vir(:, :) = 0.0d0
CALL force_env_get(force_env, potential_energy=stpcfg)
CALL plumed_gcmd_int("setStep"//CHAR(0), plumed_itimes)
CALL plumed_gcmd_doubles("setPositionsX"//CHAR(0), pos_plumed_x(:))
CALL plumed_gcmd_doubles("setPositionsY"//CHAR(0), pos_plumed_y(:))
CALL plumed_gcmd_doubles("setPositionsZ"//CHAR(0), pos_plumed_z(:))
CALL plumed_gcmd_doubles("setMasses"//CHAR(0), mass_plumed(:))
CALL plumed_gcmd_doubles("setBox"//CHAR(0), celltbt)
CALL plumed_gcmd_doubles("setVirial"//CHAR(0), plu_vir) ! PluMed Output, NOT ACTUALLY USED !
CALL plumed_gcmd_double("setEnergy"//CHAR(0), stpcfg)
CALL plumed_gcmd_doubles("setForcesX"//CHAR(0), force_plumed_x(:)) ! PluMed Output !
CALL plumed_gcmd_doubles("setForcesY"//CHAR(0), force_plumed_y(:)) ! PluMed Output !
CALL plumed_gcmd_doubles("setForcesZ"//CHAR(0), force_plumed_z(:)) ! PluMed Output !
CALL plumed_gcmd_int("prepareCalc"//CHAR(0), 0)
CALL plumed_gcmd_int("prepareDependencies"//CHAR(0), 0)
CALL plumed_gcmd_int("shareData"//CHAR(0), 0)
CALL plumed_gcmd_int("performCalc"//CHAR(0), 0)
DO i_part = 1, natom_plumed
subsys%particles%els(i_part)%f(1) = subsys%particles%els(i_part)%f(1) + force_plumed_x(i_part)
subsys%particles%els(i_part)%f(2) = subsys%particles%els(i_part)%f(2) + force_plumed_y(i_part)
subsys%particles%els(i_part)%f(3) = subsys%particles%els(i_part)%f(3) + force_plumed_z(i_part)
END DO
DEALLOCATE (pos_plumed_x, pos_plumed_y, pos_plumed_z)
DEALLOCATE (force_plumed_x, force_plumed_y, force_plumed_z)
DEALLOCATE (mass_plumed)
! Constraints ONLY of Fixed Atom type
CALL fix_atom_control(force_env)
#endif
#if !defined (__PLUMED2)
CALL cp_abort(__LOCATION__, &
"Requested to use plumed for metadynamics, but cp2k was"// &
" not compiled with plumed support.")
#endif
ELSE
IF (force_env%meta_env%langevin) THEN
IF (.NOT. PRESENT(rand)) THEN
CPABORT("Langevin on COLVAR not implemented for this MD ensemble!")
END IF
! *** Velocity Verlet for Langevin S(t)->S(t+1)
CALL metadyn_position_colvar(force_env)
! *** Forces from Vs and S(X(t+1))
CALL metadyn_forces(force_env)
! *** Velocity Verlet for Langeving *** v(t+1/2)--> v(t)
CALL metadyn_velocities_colvar(force_env, rand)
ELSE
CALL metadyn_forces(force_env, vel)
END IF
! *** Write down COVAR informations
CALL metadyn_write_colvar(force_env)
END IF
END IF
CALL timestop(handle)
END SUBROUTINE metadyn_integrator
! **************************************************************************************************
!> \brief add forces to the subsys due to the metadynamics run
!> possibly modifies the velocites (if reflective walls are applied)
!> \param force_env ...
!> \param vel ...
!> \par History
!> 04.2004 created
! **************************************************************************************************
SUBROUTINE metadyn_forces(force_env, vel)
TYPE(force_env_type), POINTER :: force_env
REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT), &
OPTIONAL :: vel
CHARACTER(len=*), PARAMETER :: routineN = 'metadyn_forces'
INTEGER :: handle, i, i_c, icolvar, ii, iwall
LOGICAL :: explicit
REAL(kind=dp) :: check_val, diff_ss, dt, ekin_w, fac_t, &
fft, norm, rval, scal, scalf, &
ss0_test, tol_ekin
TYPE(colvar_p_type), DIMENSION(:), POINTER :: colvar_p
TYPE(cp_logger_type), POINTER :: logger
TYPE(cp_subsys_type), POINTER :: subsys
TYPE(meta_env_type), POINTER :: meta_env
TYPE(metavar_type), POINTER :: cv
TYPE(particle_list_type), POINTER :: particles
TYPE(section_vals_type), POINTER :: ss0_section, vvp_section
NULLIFY (logger, meta_env)
meta_env => force_env%meta_env
IF (.NOT. ASSOCIATED(meta_env)) RETURN
CALL timeset(routineN, handle)
logger => cp_get_default_logger()
NULLIFY (colvar_p, subsys, cv, ss0_section, vvp_section)
CALL force_env_get(force_env, subsys=subsys)
dt = meta_env%dt
IF (.NOT. meta_env%restart) meta_env%n_steps = meta_env%n_steps + 1
! Initialize velocity
IF (meta_env%restart .AND. meta_env%extended_lagrange) THEN
meta_env%ekin_s = 0.0_dp
DO i_c = 1, meta_env%n_colvar
cv => meta_env%metavar(i_c)
cv%vvp = force_env%globenv%gaussian_rng_stream%next()
meta_env%ekin_s = meta_env%ekin_s + 0.5_dp*cv%mass*cv%vvp**2
END DO
ekin_w = 0.5_dp*meta_env%temp_wanted*REAL(meta_env%n_colvar, KIND=dp)
fac_t = SQRT(ekin_w/MAX(meta_env%ekin_s, 1.0E-8_dp))
DO i_c = 1, meta_env%n_colvar
cv => meta_env%metavar(i_c)
cv%vvp = cv%vvp*fac_t
END DO
meta_env%ekin_s = 0.0_dp
END IF
! *** Velocity Verlet for Langevin S(t)->S(t+1)
! compute ss and the derivative of ss with respect to the atomic positions
DO i_c = 1, meta_env%n_colvar
cv => meta_env%metavar(i_c)
icolvar = cv%icolvar
CALL colvar_eval_glob_f(icolvar, force_env)
cv%ss = subsys%colvar_p(icolvar)%colvar%ss
! Setup the periodic flag if the COLVAR is (-pi,pi] periodic
cv%periodic = (subsys%colvar_p(icolvar)%colvar%type_id == torsion_colvar_id)
! Restart for Extended Lagrangian Metadynamics
IF (meta_env%restart) THEN
! Initialize the position of the collective variable in the extended lagrange
ss0_section => section_vals_get_subs_vals(meta_env%metadyn_section, "EXT_LAGRANGE_SS0")
CALL section_vals_get(ss0_section, explicit=explicit)
IF (explicit) THEN
CALL section_vals_val_get(ss0_section, "_DEFAULT_KEYWORD_", &
i_rep_val=i_c, r_val=rval)
cv%ss0 = rval
ELSE
cv%ss0 = cv%ss
END IF
vvp_section => section_vals_get_subs_vals(meta_env%metadyn_section, "EXT_LAGRANGE_VVP")
CALL section_vals_get(vvp_section, explicit=explicit)
IF (explicit) THEN
CALL section_vals_val_get(vvp_section, "_DEFAULT_KEYWORD_", &
i_rep_val=i_c, r_val=rval)
cv%vvp = rval
END IF
END IF
!
IF (.NOT. meta_env%extended_lagrange) THEN
cv%ss0 = cv%ss
cv%vvp = 0.0_dp
END IF
END DO
! History dependent forces (evaluated at s0)
IF (meta_env%do_hills) CALL hills(meta_env)
! Apply walls to the colvars
CALL meta_walls(meta_env)
meta_env%restart = .FALSE.
IF (.NOT. meta_env%extended_lagrange) THEN
meta_env%ekin_s = 0.0_dp
meta_env%epot_s = 0.0_dp
meta_env%epot_walls = 0.0_dp
DO i_c = 1, meta_env%n_colvar
cv => meta_env%metavar(i_c)
cv%epot_s = 0.0_dp
cv%ff_s = 0.0_dp
meta_env%epot_walls = meta_env%epot_walls + cv%epot_walls
icolvar = cv%icolvar
NULLIFY (particles)
CALL cp_subsys_get(subsys, colvar_p=colvar_p, &
particles=particles)
DO ii = 1, colvar_p(icolvar)%colvar%n_atom_s
i = colvar_p(icolvar)%colvar%i_atom(ii)
fft = cv%ff_hills + cv%ff_walls
particles%els(i)%f = particles%els(i)%f + fft*colvar_p(icolvar)%colvar%dsdr(:, ii)
END DO
END DO
ELSE
meta_env%ekin_s = 0.0_dp
meta_env%epot_s = 0.0_dp
meta_env%epot_walls = 0.0_dp
DO i_c = 1, meta_env%n_colvar
cv => meta_env%metavar(i_c)
diff_ss = cv%ss - cv%ss0
IF (cv%periodic) THEN
! The difference of a periodic COLVAR is always within [-pi,pi]
diff_ss = SIGN(1.0_dp, ASIN(SIN(diff_ss)))*ACOS(COS(diff_ss))
END IF
cv%epot_s = 0.5_dp*cv%lambda*(diff_ss)**2.0_dp
cv%ff_s = cv%lambda*(diff_ss)
icolvar = cv%icolvar
! forces on the atoms
NULLIFY (particles)
CALL cp_subsys_get(subsys, colvar_p=colvar_p, &
particles=particles)
DO ii = 1, colvar_p(icolvar)%colvar%n_atom_s
i = colvar_p(icolvar)%colvar%i_atom(ii)
particles%els(i)%f = particles%els(i)%f - cv%ff_s*colvar_p(icolvar)%colvar%dsdr(:, ii)
END DO
! velocity verlet on the s0 if NOT langevin
IF (.NOT. meta_env%langevin) THEN
fft = cv%ff_s + cv%ff_hills + cv%ff_walls
cv%vvp = cv%vvp + dt*fft/cv%mass
meta_env%ekin_s = meta_env%ekin_s + 0.5_dp*cv%mass*cv%vvp**2
meta_env%epot_s = meta_env%epot_s + cv%epot_s
meta_env%epot_walls = meta_env%epot_walls + cv%epot_walls
END IF
END DO
! velocity rescaling on the s0
IF (meta_env%tempcontrol .AND. (.NOT. meta_env%langevin)) THEN
ekin_w = 0.5_dp*meta_env%temp_wanted*REAL(meta_env%n_colvar, KIND=dp)
tol_ekin = 0.5_dp*meta_env%toll_temp*REAL(meta_env%n_colvar, KIND=dp)
IF (ABS(ekin_w - meta_env%ekin_s) > tol_ekin) THEN
fac_t = SQRT(ekin_w/MAX(meta_env%ekin_s, 1.0E-8_dp))
DO i_c = 1, meta_env%n_colvar
cv => meta_env%metavar(i_c)
cv%vvp = cv%vvp*fac_t
END DO
meta_env%ekin_s = ekin_w
END IF
END IF
! Reflective Wall only for s0
DO i_c = 1, meta_env%n_colvar
cv => meta_env%metavar(i_c)
IF (cv%do_wall) THEN
DO iwall = 1, SIZE(cv%walls)
SELECT CASE (cv%walls(iwall)%id_type)
CASE (do_wall_reflective)
ss0_test = cv%ss0 + dt*cv%vvp
IF (cv%periodic) THEN
! A periodic COLVAR is always within [-pi,pi]
ss0_test = SIGN(1.0_dp, ASIN(SIN(ss0_test)))*ACOS(COS(ss0_test))
END IF
SELECT CASE (cv%walls(iwall)%id_direction)
CASE (do_wall_p)
IF ((ss0_test > cv%walls(iwall)%pos) .AND. (cv%vvp > 0)) cv%vvp = -cv%vvp
CASE (do_wall_m)
IF ((ss0_test < cv%walls(iwall)%pos) .AND. (cv%vvp < 0)) cv%vvp = -cv%vvp
END SELECT
END SELECT
END DO
END IF
END DO
! Update of ss0 if NOT langevin
IF (.NOT. meta_env%langevin) THEN
DO i_c = 1, meta_env%n_colvar
cv => meta_env%metavar(i_c)
cv%ss0 = cv%ss0 + dt*cv%vvp
IF (cv%periodic) THEN
! A periodic COLVAR is always within [-pi,pi]
cv%ss0 = SIGN(1.0_dp, ASIN(SIN(cv%ss0)))*ACOS(COS(cv%ss0))
END IF
END DO
END IF
END IF
! Constraints ONLY of Fixed Atom type
CALL fix_atom_control(force_env)
! Reflective Wall only for ss
DO i_c = 1, meta_env%n_colvar
cv => meta_env%metavar(i_c)
IF (cv%do_wall) THEN
DO iwall = 1, SIZE(cv%walls)
SELECT CASE (cv%walls(iwall)%id_type)
CASE (do_wall_reflective)
SELECT CASE (cv%walls(iwall)%id_direction)
CASE (do_wall_p)
IF (cv%ss < cv%walls(iwall)%pos) CYCLE
check_val = -1.0_dp
CASE (do_wall_m)
IF (cv%ss > cv%walls(iwall)%pos) CYCLE
check_val = 1.0_dp
END SELECT
NULLIFY (particles)
icolvar = cv%icolvar
CALL cp_subsys_get(subsys, colvar_p=colvar_p, particles=particles)
scal = 0.0_dp
scalf = 0.0_dp
norm = 0.0_dp
DO ii = 1, colvar_p(icolvar)%colvar%n_atom_s
i = colvar_p(icolvar)%colvar%i_atom(ii)
IF (PRESENT(vel)) THEN
scal = scal + vel(1, i)*colvar_p(icolvar)%colvar%dsdr(1, ii)
scal = scal + vel(2, i)*colvar_p(icolvar)%colvar%dsdr(2, ii)
scal = scal + vel(3, i)*colvar_p(icolvar)%colvar%dsdr(3, ii)
ELSE
scal = scal + particles%els(i)%v(1)*colvar_p(icolvar)%colvar%dsdr(1, ii)
scal = scal + particles%els(i)%v(2)*colvar_p(icolvar)%colvar%dsdr(2, ii)
scal = scal + particles%els(i)%v(3)*colvar_p(icolvar)%colvar%dsdr(3, ii)
END IF
scalf = scalf + particles%els(i)%f(1)*colvar_p(icolvar)%colvar%dsdr(1, ii)
scalf = scalf + particles%els(i)%f(2)*colvar_p(icolvar)%colvar%dsdr(2, ii)
scalf = scalf + particles%els(i)%f(3)*colvar_p(icolvar)%colvar%dsdr(3, ii)
norm = norm + colvar_p(icolvar)%colvar%dsdr(1, ii)**2
norm = norm + colvar_p(icolvar)%colvar%dsdr(2, ii)**2
norm = norm + colvar_p(icolvar)%colvar%dsdr(3, ii)**2
END DO
IF (norm /= 0.0_dp) scal = scal/norm
IF (norm /= 0.0_dp) scalf = scalf/norm
IF (scal*check_val > 0.0_dp) CYCLE
DO ii = 1, colvar_p(icolvar)%colvar%n_atom_s
i = colvar_p(icolvar)%colvar%i_atom(ii)
IF (PRESENT(vel)) THEN
vel(:, i) = vel(:, i) - 2.0_dp*colvar_p(icolvar)%colvar%dsdr(:, ii)*scal
ELSE
particles%els(i)%v(:) = particles%els(i)%v(:) - 2.0_dp*colvar_p(icolvar)%colvar%dsdr(:, ii)*scal
END IF
! Nullify forces along the colvar (this avoids the weird behaviors of the reflective wall)
particles%els(i)%f(:) = particles%els(i)%f(:) - colvar_p(icolvar)%colvar%dsdr(:, ii)*scalf
END DO
END SELECT
END DO
END IF
END DO
CALL timestop(handle)
END SUBROUTINE metadyn_forces
! **************************************************************************************************
!> \brief Evolves velocities COLVAR according to
!> Vanden-Eijnden Ciccotti C.Phys.Letter 429 (2006) 310-316
!> \param force_env ...
!> \param rand ...
!> \date 01.2009
!> \author Fabio Sterpone and Teodoro Laino
! **************************************************************************************************
SUBROUTINE metadyn_velocities_colvar(force_env, rand)
TYPE(force_env_type), POINTER :: force_env
REAL(KIND=dp), DIMENSION(:), INTENT(INOUT), &
OPTIONAL :: rand
CHARACTER(len=*), PARAMETER :: routineN = 'metadyn_velocities_colvar'
INTEGER :: handle, i_c
REAL(kind=dp) :: diff_ss, dt, fft, sigma
TYPE(cp_logger_type), POINTER :: logger
TYPE(meta_env_type), POINTER :: meta_env
TYPE(metavar_type), POINTER :: cv
NULLIFY (logger, meta_env, cv)
meta_env => force_env%meta_env
IF (.NOT. ASSOCIATED(meta_env)) RETURN
CALL timeset(routineN, handle)
logger => cp_get_default_logger()
! Add citation
IF (meta_env%langevin) CALL cite_reference(VandenCic2006)
dt = meta_env%dt
! History dependent forces (evaluated at s0)
IF (meta_env%do_hills) CALL hills(meta_env)
! Evolve Velocities
meta_env%ekin_s = 0.0_dp
meta_env%epot_walls = 0.0_dp
DO i_c = 1, meta_env%n_colvar
cv => meta_env%metavar(i_c)
diff_ss = cv%ss - cv%ss0
IF (cv%periodic) THEN
! The difference of a periodic COLVAR is always within [-pi,pi]
diff_ss = SIGN(1.0_dp, ASIN(SIN(diff_ss)))*ACOS(COS(diff_ss))
END IF
cv%epot_s = 0.5_dp*cv%lambda*(diff_ss)**2.0_dp
cv%ff_s = cv%lambda*(diff_ss)
fft = cv%ff_s + cv%ff_hills
sigma = SQRT((meta_env%temp_wanted*kelvin)*2.0_dp*(boltzmann/joule)*cv%gamma/cv%mass)
cv%vvp = cv%vvp + 0.5_dp*dt*fft/cv%mass - 0.5_dp*dt*cv%gamma*cv%vvp + &
0.5_dp*SQRT(dt)*sigma*rand(i_c)
meta_env%ekin_s = meta_env%ekin_s + 0.5_dp*cv%mass*cv%vvp**2
meta_env%epot_walls = meta_env%epot_walls + cv%epot_walls
END DO
CALL timestop(handle)
END SUBROUTINE metadyn_velocities_colvar
! **************************************************************************************************
!> \brief Evolves COLVAR position
!> Vanden-Eijnden Ciccotti C.Phys.Letter 429 (2006) 310-316
!> \param force_env ...
!> \date 01.2009
!> \author Fabio Sterpone and Teodoro Laino
! **************************************************************************************************
SUBROUTINE metadyn_position_colvar(force_env)
TYPE(force_env_type), POINTER :: force_env
CHARACTER(len=*), PARAMETER :: routineN = 'metadyn_position_colvar'
INTEGER :: handle, i_c
REAL(kind=dp) :: dt
TYPE(cp_logger_type), POINTER :: logger
TYPE(meta_env_type), POINTER :: meta_env
TYPE(metavar_type), POINTER :: cv
NULLIFY (logger, meta_env, cv)
meta_env => force_env%meta_env
IF (.NOT. ASSOCIATED(meta_env)) RETURN
CALL timeset(routineN, handle)
logger => cp_get_default_logger()
! Add citation
IF (meta_env%langevin) CALL cite_reference(VandenCic2006)
dt = meta_env%dt
! Update of ss0
DO i_c = 1, meta_env%n_colvar
cv => meta_env%metavar(i_c)
cv%ss0 = cv%ss0 + dt*cv%vvp
IF (cv%periodic) THEN
! A periodic COLVAR is always within [-pi,pi]
cv%ss0 = SIGN(1.0_dp, ASIN(SIN(cv%ss0)))*ACOS(COS(cv%ss0))
END IF
END DO
CALL timestop(handle)
END SUBROUTINE metadyn_position_colvar
! **************************************************************************************************
!> \brief Write down COLVAR information evolved according to
!> Vanden-Eijnden Ciccotti C.Phys.Letter 429 (2006) 310-316
!> \param force_env ...
!> \date 01.2009
!> \author Fabio Sterpone and Teodoro Laino
! **************************************************************************************************
SUBROUTINE metadyn_write_colvar(force_env)
TYPE(force_env_type), POINTER :: force_env
CHARACTER(len=*), PARAMETER :: routineN = 'metadyn_write_colvar'
INTEGER :: handle, i, i_c, iw
REAL(KIND=dp) :: diff_ss, temp
TYPE(cp_logger_type), POINTER :: logger
TYPE(meta_env_type), POINTER :: meta_env
TYPE(metavar_type), POINTER :: cv
NULLIFY (logger, meta_env, cv)
meta_env => force_env%meta_env
IF (.NOT. ASSOCIATED(meta_env)) RETURN
CALL timeset(routineN, handle)
logger => cp_get_default_logger()
! If Langevin we need to recompute few quantities
! This does not apply to the standard lagrangian scheme since it is
! implemented with a plain Newton integration scheme.. while Langevin
! follows the correct Verlet integration.. This will have to be made
! uniform in the future (Teodoro Laino - 01.2009)
IF (meta_env%langevin) THEN
meta_env%ekin_s = 0.0_dp
meta_env%epot_s = 0.0_dp
DO i_c = 1, meta_env%n_colvar
cv => meta_env%metavar(i_c)
diff_ss = cv%ss - cv%ss0
IF (cv%periodic) THEN
! The difference of a periodic COLVAR is always within [-pi,pi]
diff_ss = SIGN(1.0_dp, ASIN(SIN(diff_ss)))*ACOS(COS(diff_ss))
END IF
cv%epot_s = 0.5_dp*cv%lambda*(diff_ss)**2.0_dp
cv%ff_s = cv%lambda*(diff_ss)
meta_env%epot_s = meta_env%epot_s + cv%epot_s
meta_env%ekin_s = meta_env%ekin_s + 0.5_dp*cv%mass*cv%vvp**2
END DO
END IF
! write COLVAR file
iw = cp_print_key_unit_nr(logger, meta_env%metadyn_section, &
"PRINT%COLVAR", extension=".metadynLog")
IF (iw > 0) THEN
IF (meta_env%extended_lagrange) THEN
WRITE (iw, '(f16.8,70f15.8)') meta_env%time*femtoseconds, &
(meta_env%metavar(i)%ss0, i=1, meta_env%n_colvar), &
(meta_env%metavar(i)%ss, i=1, meta_env%n_colvar), &
(meta_env%metavar(i)%ff_s, i=1, meta_env%n_colvar), &
(meta_env%metavar(i)%ff_hills, i=1, meta_env%n_colvar), &
(meta_env%metavar(i)%ff_walls, i=1, meta_env%n_colvar), &
(meta_env%metavar(i)%vvp, i=1, meta_env%n_colvar), &
meta_env%epot_s, &
meta_env%hills_env%energy, &
meta_env%epot_walls, &
(meta_env%ekin_s)*2.0_dp/(REAL(meta_env%n_colvar, KIND=dp))*kelvin
ELSE
WRITE (iw, '(f16.8,40f13.5)') meta_env%time*femtoseconds, &
(meta_env%metavar(i)%ss0, i=1, meta_env%n_colvar), &
(meta_env%metavar(i)%ff_hills, i=1, meta_env%n_colvar), &
(meta_env%metavar(i)%ff_walls, i=1, meta_env%n_colvar), &
meta_env%hills_env%energy, &
meta_env%epot_walls
END IF
END IF
CALL cp_print_key_finished_output(iw, logger, meta_env%metadyn_section, &
"PRINT%COLVAR")
! Temperature for COLVAR
IF (meta_env%extended_lagrange) THEN
temp = meta_env%ekin_s*2.0_dp/(REAL(meta_env%n_colvar, KIND=dp))*kelvin
meta_env%avg_temp = (meta_env%avg_temp*REAL(meta_env%n_steps, KIND=dp) + &
temp)/REAL(meta_env%n_steps + 1, KIND=dp)
iw = cp_print_key_unit_nr(logger, meta_env%metadyn_section, &
"PRINT%TEMPERATURE_COLVAR", extension=".metadynLog")
IF (iw > 0) THEN
WRITE (iw, '(T2,79("-"))')
WRITE (iw, '( A,T51,f10.2,T71,f10.2)') ' COLVARS INSTANTANEOUS/AVERAGE TEMPERATURE ', &
temp, meta_env%avg_temp
WRITE (iw, '(T2,79("-"))')
END IF
CALL cp_print_key_finished_output(iw, logger, meta_env%metadyn_section, &
"PRINT%TEMPERATURE_COLVAR")
END IF
CALL timestop(handle)
END SUBROUTINE metadyn_write_colvar
! **************************************************************************************************
!> \brief Major driver for adding hills and computing forces due to the history
!> dependent term
!> \param meta_env ...
!> \par History
!> 04.2004 created
!> 10.2008 Teodoro Laino [tlaino] - University of Zurich
!> Major rewriting and addition of multiple walkers
! **************************************************************************************************
SUBROUTINE hills(meta_env)
TYPE(meta_env_type), POINTER :: meta_env
CHARACTER(len=*), PARAMETER :: routineN = 'hills'
INTEGER :: handle, i, i_hills, ih, intermeta_steps, &
iter_nr, iw, n_colvar, n_hills_start, &
n_step
LOGICAL :: force_gauss
REAL(KIND=dp) :: cut_func, dfunc, diff, dp2, frac, &
slow_growth, V_now_here, V_to_fes, &
wtww, ww
REAL(KIND=dp), DIMENSION(:), POINTER :: ddp, denf, diff_ss, local_last_hills, &
numf
TYPE(cp_logger_type), POINTER :: logger
TYPE(hills_env_type), POINTER :: hills_env
TYPE(metavar_type), DIMENSION(:), POINTER :: colvars
TYPE(multiple_walkers_type), POINTER :: multiple_walkers
CALL timeset(routineN, handle)
NULLIFY (hills_env, multiple_walkers, logger, colvars, ddp, local_last_hills)
hills_env => meta_env%hills_env
logger => cp_get_default_logger()
colvars => meta_env%metavar
n_colvar = meta_env%n_colvar
n_step = meta_env%n_steps
! Create a temporary logger level specific for metadynamics
CALL cp_add_iter_level(logger%iter_info, "METADYNAMICS")
CALL get_meta_iter_level(meta_env, iter_nr)
CALL cp_iterate(logger%iter_info, last=.FALSE., iter_nr=iter_nr)
! Set-up restart if any
IF (meta_env%hills_env%restart) THEN
meta_env%hills_env%restart = .FALSE.
IF (meta_env%well_tempered) THEN
CALL restart_hills(hills_env%ss_history, hills_env%delta_s_history, hills_env%ww_history, &
hills_env%ww, hills_env%n_hills, n_colvar, colvars, meta_env%metadyn_section, &
invdt_history=hills_env%invdt_history)
ELSE
CALL restart_hills(hills_env%ss_history, hills_env%delta_s_history, hills_env%ww_history, &
hills_env%ww, hills_env%n_hills, n_colvar, colvars, meta_env%metadyn_section)
END IF
END IF
! Proceed with normal calculation
intermeta_steps = n_step - hills_env%old_hill_step
force_gauss = .FALSE.
IF ((hills_env%min_disp > 0.0_dp) .AND. (hills_env%old_hill_number > 0) .AND. &
(intermeta_steps >= hills_env%min_nt_hills)) THEN
ALLOCATE (ddp(meta_env%n_colvar))
ALLOCATE (local_last_hills(meta_env%n_colvar))
local_last_hills(1:n_colvar) = hills_env%ss_history(1:n_colvar, hills_env%old_hill_number)
!RG Calculate the displacement
dp2 = 0.0_dp
DO i = 1, n_colvar
ddp(i) = colvars(i)%ss0 - local_last_hills(i)
IF (colvars(i)%periodic) THEN
! The difference of a periodic COLVAR is always within [-pi,pi]
ddp(i) = SIGN(1.0_dp, ASIN(SIN(ddp(i))))*ACOS(COS(ddp(i)))
END IF
dp2 = dp2 + ddp(i)**2
END DO
dp2 = SQRT(dp2)
IF (dp2 > hills_env%min_disp) THEN
force_gauss = .TRUE.
iw = cp_print_key_unit_nr(logger, meta_env%metadyn_section, &
"PRINT%PROGRAM_RUN_INFO", extension=".metadynLog")
IF (iw > 0) THEN
WRITE (UNIT=iw, FMT="(A,F0.6,A,F0.6)") &
" METADYN| Threshold value for COLVAR displacement exceeded: ", &
dp2, " > ", hills_env%min_disp
END IF
CALL cp_print_key_finished_output(iw, logger, meta_env%metadyn_section, &
"PRINT%PROGRAM_RUN_INFO")
END IF
DEALLOCATE (ddp)
DEALLOCATE (local_last_hills)
END IF
!RG keep into account adaptive hills
IF (((MODULO(intermeta_steps, hills_env%nt_hills) == 0) .OR. force_gauss) &
.AND. (.NOT. meta_env%restart) .AND. (hills_env%nt_hills > 0)) THEN
IF (meta_env%do_multiple_walkers) multiple_walkers => meta_env%multiple_walkers
n_hills_start = hills_env%n_hills
! Add the hill corresponding to this location
IF (meta_env%well_tempered) THEN
! Well-Tempered scaling of hills height
V_now_here = 0._dp
DO ih = 1, hills_env%n_hills
dp2 = 0._dp
DO i = 1, n_colvar
diff = colvars(i)%ss0 - hills_env%ss_history(i, ih)
IF (colvars(i)%periodic) THEN
! The difference of a periodic COLVAR is always within [-pi,pi]
diff = SIGN(1.0_dp, ASIN(SIN(diff)))*ACOS(COS(diff))
END IF
diff = (diff)/hills_env%delta_s_history(i, ih)
dp2 = dp2 + diff**2
END DO
V_to_fes = 1.0_dp + meta_env%wttemperature*hills_env%invdt_history(ih)
V_now_here = V_now_here + hills_env%ww_history(ih)/V_to_fes*EXP(-0.5_dp*dp2)
END DO
wtww = hills_env%ww*EXP(-V_now_here*meta_env%invdt)
ww = wtww*(1.0_dp + meta_env%wttemperature*meta_env%invdt)
CALL add_hill_single(hills_env, colvars, ww, hills_env%n_hills, n_colvar, meta_env%invdt)
ELSE
CALL add_hill_single(hills_env, colvars, hills_env%ww, hills_env%n_hills, n_colvar)
END IF
! Update local n_hills counter
IF (meta_env%do_multiple_walkers) multiple_walkers%n_hills_local = multiple_walkers%n_hills_local + 1
hills_env%old_hill_number = hills_env%n_hills
hills_env%old_hill_step = n_step
! Update iteration level for printing
CALL get_meta_iter_level(meta_env, iter_nr)
CALL cp_iterate(logger%iter_info, last=.FALSE., iter_nr=iter_nr)
! Print just program_run_info
iw = cp_print_key_unit_nr(logger, meta_env%metadyn_section, &
"PRINT%PROGRAM_RUN_INFO", extension=".metadynLog")
IF (iw > 0) THEN
IF (meta_env%do_multiple_walkers) THEN
WRITE (iw, '(/,1X,"METADYN|",A,I0,A,I0,A,/)') &
' Global/Local Hills number (', hills_env%n_hills, '/', multiple_walkers%n_hills_local, &
') added.'
ELSE
WRITE (iw, '(/,1X,"METADYN|",A,I0,A,/)') ' Hills number (', hills_env%n_hills, ') added.'
END IF
END IF
CALL cp_print_key_finished_output(iw, logger, meta_env%metadyn_section, &
"PRINT%PROGRAM_RUN_INFO")
! Handle Multiple Walkers
IF (meta_env%do_multiple_walkers) THEN
! Print Local Hills file if requested
iw = cp_print_key_unit_nr(logger, meta_env%metadyn_section, &
"PRINT%HILLS", middle_name="LOCAL", extension=".metadynLog")
IF (iw > 0) THEN
WRITE (iw, '(f12.1,30f13.5)') meta_env%time*femtoseconds, &
(hills_env%ss_history(ih, hills_env%n_hills), ih=1, n_colvar), &
(hills_env%delta_s_history(ih, hills_env%n_hills), ih=1, n_colvar), &
hills_env%ww_history(hills_env%n_hills)
END IF
CALL cp_print_key_finished_output(iw, logger, meta_env%metadyn_section, &
"PRINT%HILLS")
! Check the communication buffer of the other walkers
CALL synchronize_multiple_walkers(multiple_walkers, hills_env, colvars, &
n_colvar, meta_env%metadyn_section)
END IF
! Print Hills file if requested (for multiple walkers this file includes
! the Hills coming from all the walkers).
iw = cp_print_key_unit_nr(logger, meta_env%metadyn_section, &
"PRINT%HILLS", extension=".metadynLog")
IF (iw > 0) THEN
DO i_hills = n_hills_start + 1, hills_env%n_hills
WRITE (iw, '(f12.1,30f13.5)') meta_env%time*femtoseconds, &
(hills_env%ss_history(ih, i_hills), ih=1, n_colvar), &
(hills_env%delta_s_history(ih, i_hills), ih=1, n_colvar), &
hills_env%ww_history(i_hills)
END DO
END IF
CALL cp_print_key_finished_output(iw, logger, meta_env%metadyn_section, &
"PRINT%HILLS")
END IF
! Computes forces due to the hills: history dependent term