-
Notifications
You must be signed in to change notification settings - Fork 1
/
fist_pol_scf.F
729 lines (673 loc) · 37.5 KB
/
fist_pol_scf.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
!--------------------------------------------------------------------------------------------------!
! 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 !
!--------------------------------------------------------------------------------------------------!
! **************************************************************************************************
!> \par History
!> CJM APRIL-30-2009: only uses fist_env
!> Teodoro Laino [tlaino] - 05.2009 : Generalization to different Ewald
!> methods (initial framework)
!> \author CJM
! **************************************************************************************************
MODULE fist_pol_scf
USE atomic_kind_types, ONLY: atomic_kind_type,&
get_atomic_kind
USE cell_types, ONLY: cell_type
USE cp_log_handling, ONLY: cp_get_default_logger,&
cp_logger_type
USE cp_output_handling, ONLY: cp_print_key_finished_output,&
cp_print_key_unit_nr
USE distribution_1d_types, ONLY: distribution_1d_type
USE ewald_environment_types, ONLY: ewald_env_get,&
ewald_environment_type
USE ewald_pw_types, ONLY: ewald_pw_type
USE ewalds_multipole, ONLY: ewald_multipole_evaluate
USE fist_energy_types, ONLY: fist_energy_type
USE fist_nonbond_env_types, ONLY: fist_nonbond_env_type
USE input_constants, ONLY: do_fist_pol_cg,&
do_fist_pol_sc
USE input_section_types, ONLY: section_vals_type
USE kinds, ONLY: dp
USE multipole_types, ONLY: multipole_type
USE particle_types, ONLY: particle_type
USE pw_poisson_types, ONLY: do_ewald_ewald,&
do_ewald_pme,&
do_ewald_spme
#include "./base/base_uses.f90"
IMPLICIT NONE
PRIVATE
LOGICAL, PRIVATE :: debug_this_module = .FALSE.
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'fist_pol_scf'
PUBLIC :: fist_pol_evaluate
CONTAINS
! **************************************************************************************************
!> \brief Main driver for evaluating energy/forces in a polarizable forcefield
!> \param atomic_kind_set ...
!> \param multipoles ...
!> \param ewald_env ...
!> \param ewald_pw ...
!> \param fist_nonbond_env ...
!> \param cell ...
!> \param particle_set ...
!> \param local_particles ...
!> \param thermo ...
!> \param vg_coulomb ...
!> \param pot_nonbond ...
!> \param f_nonbond ...
!> \param fg_coulomb ...
!> \param use_virial ...
!> \param pv_g ...
!> \param pv_nonbond ...
!> \param mm_section ...
!> \param do_ipol ...
!> \author [email protected] (2010-03-01)
! **************************************************************************************************
SUBROUTINE fist_pol_evaluate(atomic_kind_set, multipoles, ewald_env, &
ewald_pw, fist_nonbond_env, cell, particle_set, local_particles, &
thermo, vg_coulomb, pot_nonbond, f_nonbond, fg_coulomb, use_virial, &
pv_g, pv_nonbond, mm_section, do_ipol)
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(multipole_type), POINTER :: multipoles
TYPE(ewald_environment_type), POINTER :: ewald_env
TYPE(ewald_pw_type), POINTER :: ewald_pw
TYPE(fist_nonbond_env_type), POINTER :: fist_nonbond_env
TYPE(cell_type), POINTER :: cell
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(distribution_1d_type), POINTER :: local_particles
TYPE(fist_energy_type), POINTER :: thermo
REAL(KIND=dp) :: vg_coulomb, pot_nonbond
REAL(KIND=dp), DIMENSION(:, :) :: f_nonbond, fg_coulomb
LOGICAL, INTENT(IN) :: use_virial
REAL(KIND=dp), DIMENSION(3, 3) :: pv_g, pv_nonbond
TYPE(section_vals_type), POINTER :: mm_section
INTEGER :: do_ipol
SELECT CASE (do_ipol)
CASE (do_fist_pol_sc)
CALL fist_pol_evaluate_sc(atomic_kind_set, multipoles, ewald_env, &
ewald_pw, fist_nonbond_env, cell, particle_set, local_particles, &
thermo, vg_coulomb, pot_nonbond, f_nonbond, fg_coulomb, use_virial, &
pv_g, pv_nonbond, mm_section)
CASE (do_fist_pol_cg)
CALL fist_pol_evaluate_cg(atomic_kind_set, multipoles, ewald_env, &
ewald_pw, fist_nonbond_env, cell, particle_set, local_particles, &
thermo, vg_coulomb, pot_nonbond, f_nonbond, fg_coulomb, use_virial, &
pv_g, pv_nonbond, mm_section)
END SELECT
END SUBROUTINE fist_pol_evaluate
! **************************************************************************************************
!> \brief Self-consistent solver for a polarizable force-field
!> \param atomic_kind_set ...
!> \param multipoles ...
!> \param ewald_env ...
!> \param ewald_pw ...
!> \param fist_nonbond_env ...
!> \param cell ...
!> \param particle_set ...
!> \param local_particles ...
!> \param thermo ...
!> \param vg_coulomb ...
!> \param pot_nonbond ...
!> \param f_nonbond ...
!> \param fg_coulomb ...
!> \param use_virial ...
!> \param pv_g ...
!> \param pv_nonbond ...
!> \param mm_section ...
!> \author [email protected] (2010-03-01)
!> \note
!> Method: Given an initial guess of the induced dipoles, the electrostatic
!> field is computed at each dipole. Then new induced dipoles are computed
!> following p = alpha x E. This is repeated until a convergence criterion is
!> met. The convergence is measured as the RSMD of the derivatives of the
!> electrostatic energy (including dipole self-energy) towards the components
!> of the dipoles.
! **************************************************************************************************
SUBROUTINE fist_pol_evaluate_sc(atomic_kind_set, multipoles, ewald_env, ewald_pw, &
fist_nonbond_env, cell, particle_set, local_particles, thermo, vg_coulomb, &
pot_nonbond, f_nonbond, fg_coulomb, use_virial, pv_g, pv_nonbond, mm_section)
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(multipole_type), POINTER :: multipoles
TYPE(ewald_environment_type), POINTER :: ewald_env
TYPE(ewald_pw_type), POINTER :: ewald_pw
TYPE(fist_nonbond_env_type), POINTER :: fist_nonbond_env
TYPE(cell_type), POINTER :: cell
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(distribution_1d_type), POINTER :: local_particles
TYPE(fist_energy_type), POINTER :: thermo
REAL(KIND=dp) :: vg_coulomb, pot_nonbond
REAL(KIND=dp), DIMENSION(:, :) :: f_nonbond, fg_coulomb
LOGICAL, INTENT(IN) :: use_virial
REAL(KIND=dp), DIMENSION(3, 3) :: pv_g, pv_nonbond
TYPE(section_vals_type), POINTER :: mm_section
CHARACTER(len=*), PARAMETER :: routineN = 'fist_pol_evaluate_sc'
INTEGER :: ewald_type, handle, i, iatom, ii, ikind, &
iter, iw, iw2, j, max_ipol_iter, &
natom_of_kind, natoms, nkind, ntot
INTEGER, DIMENSION(:), POINTER :: atom_list
LOGICAL :: iwarn
REAL(KIND=dp) :: apol, cpol, eps_pol, pot_nonbond_local, &
rmsd, tmp_trace
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: efield1, efield2
TYPE(atomic_kind_type), POINTER :: atomic_kind
TYPE(cp_logger_type), POINTER :: logger
CALL timeset(routineN, handle)
NULLIFY (logger, atomic_kind)
logger => cp_get_default_logger()
iw = cp_print_key_unit_nr(logger, mm_section, "PRINT%ITER_INFO", &
extension=".mmLog")
iw2 = cp_print_key_unit_nr(logger, mm_section, "PRINT%EWALD_INFO", &
extension=".mmLog")
CALL ewald_env_get(ewald_env, max_ipol_iter=max_ipol_iter, eps_pol=eps_pol, &
ewald_type=ewald_type)
natoms = SIZE(particle_set)
ALLOCATE (efield1(3, natoms))
ALLOCATE (efield2(9, natoms))
nkind = SIZE(atomic_kind_set)
IF (iw > 0) THEN
WRITE (iw, FMT='(/,T2,"POL_SCF|" ,"Method: self-consistent")')
WRITE (iw, FMT='(T2,"POL_SCF| "," Iteration",7X,"Conv.",T49,"Electrostatic & Induction Energy")')
END IF
pol_scf: DO iter = 1, max_ipol_iter
! Evaluate the electrostatic with Ewald schemes
CALL eval_pol_ewald(ewald_type, ewald_env, ewald_pw, fist_nonbond_env, cell, &
particle_set, local_particles, vg_coulomb, pot_nonbond_local, thermo, &
multipoles, do_correction_bonded=.TRUE., do_forces=.FALSE., &
do_stress=.FALSE., do_efield=.TRUE., iw2=iw2, do_debug=.FALSE., &
atomic_kind_set=atomic_kind_set, mm_section=mm_section, &
efield1=efield1, efield2=efield2)
CALL logger%para_env%sum(pot_nonbond_local)
! compute the new dipoles, qudrupoles, and check for convergence
ntot = 0
rmsd = 0.0_dp
thermo%e_induction = 0.0_dp
DO ikind = 1, nkind
atomic_kind => atomic_kind_set(ikind)
CALL get_atomic_kind(atomic_kind, apol=apol, cpol=cpol, natom=natom_of_kind, atom_list=atom_list)
! ignore atoms with dipole and quadrupole polarizability zero
IF (apol == 0 .AND. cpol == 0) CYCLE
! increment counter correctly
IF (apol /= 0) ntot = ntot + natom_of_kind
IF (cpol /= 0) ntot = ntot + natom_of_kind
DO iatom = 1, natom_of_kind
ii = atom_list(iatom)
IF (apol /= 0) THEN
DO i = 1, 3
! the rmsd of the derivatives of the energy towards the
! components of the atomic dipole moments
rmsd = rmsd + (multipoles%dipoles(i, ii)/apol - efield1(i, ii))**2
END DO
END IF
IF (cpol /= 0) THEN
rmsd = rmsd + (multipoles%quadrupoles(1, 1, ii)/cpol - efield2(1, ii))**2
rmsd = rmsd + (multipoles%quadrupoles(2, 1, ii)/cpol - efield2(2, ii))**2
rmsd = rmsd + (multipoles%quadrupoles(3, 1, ii)/cpol - efield2(3, ii))**2
rmsd = rmsd + (multipoles%quadrupoles(1, 2, ii)/cpol - efield2(4, ii))**2
rmsd = rmsd + (multipoles%quadrupoles(2, 2, ii)/cpol - efield2(5, ii))**2
rmsd = rmsd + (multipoles%quadrupoles(3, 2, ii)/cpol - efield2(6, ii))**2
rmsd = rmsd + (multipoles%quadrupoles(1, 3, ii)/cpol - efield2(7, ii))**2
rmsd = rmsd + (multipoles%quadrupoles(2, 3, ii)/cpol - efield2(8, ii))**2
rmsd = rmsd + (multipoles%quadrupoles(3, 3, ii)/cpol - efield2(9, ii))**2
END IF
! compute dipole
multipoles%dipoles(:, ii) = apol*efield1(:, ii)
! compute quadrupole
IF (cpol /= 0) THEN
multipoles%quadrupoles(1, 1, ii) = cpol*efield2(1, ii)
multipoles%quadrupoles(2, 1, ii) = cpol*efield2(2, ii)
multipoles%quadrupoles(3, 1, ii) = cpol*efield2(3, ii)
multipoles%quadrupoles(1, 2, ii) = cpol*efield2(4, ii)
multipoles%quadrupoles(2, 2, ii) = cpol*efield2(5, ii)
multipoles%quadrupoles(3, 2, ii) = cpol*efield2(6, ii)
multipoles%quadrupoles(1, 3, ii) = cpol*efield2(7, ii)
multipoles%quadrupoles(2, 3, ii) = cpol*efield2(8, ii)
multipoles%quadrupoles(3, 3, ii) = cpol*efield2(9, ii)
END IF
! Compute the new induction term while we are here
IF (apol /= 0) THEN
thermo%e_induction = thermo%e_induction + &
DOT_PRODUCT(multipoles%dipoles(:, ii), &
multipoles%dipoles(:, ii))/apol/2.0_dp
END IF
IF (cpol /= 0) THEN
tmp_trace = 0._dp
DO i = 1, 3
DO j = 1, 3
tmp_trace = tmp_trace + &
multipoles%quadrupoles(i, j, ii)*multipoles%quadrupoles(i, j, ii)
END DO
END DO
thermo%e_induction = thermo%e_induction + tmp_trace/cpol/6.0_dp
END IF
END DO
END DO
rmsd = SQRT(rmsd/REAL(ntot, KIND=dp))
IF (iw > 0) THEN
! print the energy that is minimized (this is electrostatic + induction)
WRITE (iw, FMT='(T5,"POL_SCF|",5X,I5,5X,E12.6,T61,F20.10)') iter, &
rmsd, vg_coulomb + pot_nonbond_local + thermo%e_induction
END IF
IF (rmsd <= eps_pol) THEN
IF (iw > 0) WRITE (iw, FMT='(T5,"POL_SCF|",1X,"Self-consistent Polarization achieved.")')
EXIT pol_scf
END IF
iwarn = ((rmsd > eps_pol) .AND. (iter == max_ipol_iter))
IF (iwarn .AND. iw > 0) WRITE (iw, FMT='(T5,"POL_SCF|",1X,"Self-consistent Polarization not converged!")')
CPWARN_IF(iwarn, "Self-consistent Polarization not converged!")
END DO pol_scf
! Now evaluate after convergence to obtain forces and converged energies
CALL eval_pol_ewald(ewald_type, ewald_env, ewald_pw, fist_nonbond_env, cell, &
particle_set, local_particles, vg_coulomb, pot_nonbond_local, thermo, &
multipoles, do_correction_bonded=.TRUE., do_forces=.TRUE., &
do_stress=use_virial, do_efield=.FALSE., iw2=iw2, do_debug=.FALSE., &
atomic_kind_set=atomic_kind_set, mm_section=mm_section, &
forces_local=fg_coulomb, forces_glob=f_nonbond, &
pv_local=pv_g, pv_glob=pv_nonbond)
pot_nonbond = pot_nonbond + pot_nonbond_local
CALL logger%para_env%sum(pot_nonbond_local)
IF (iw > 0) THEN
! print the energy that is minimized (this is electrostatic + induction)
WRITE (iw, FMT='(T5,"POL_SCF|",5X,"Final",T61,F20.10,/)') &
vg_coulomb + pot_nonbond_local + thermo%e_induction
END IF
! Deallocate working arrays
DEALLOCATE (efield1)
DEALLOCATE (efield2)
CALL cp_print_key_finished_output(iw2, logger, mm_section, &
"PRINT%EWALD_INFO")
CALL cp_print_key_finished_output(iw, logger, mm_section, &
"PRINT%ITER_INFO")
CALL timestop(handle)
END SUBROUTINE fist_pol_evaluate_sc
! **************************************************************************************************
!> \brief Conjugate-gradient solver for a polarizable force-field
!> \param atomic_kind_set ...
!> \param multipoles ...
!> \param ewald_env ...
!> \param ewald_pw ...
!> \param fist_nonbond_env ...
!> \param cell ...
!> \param particle_set ...
!> \param local_particles ...
!> \param thermo ...
!> \param vg_coulomb ...
!> \param pot_nonbond ...
!> \param f_nonbond ...
!> \param fg_coulomb ...
!> \param use_virial ...
!> \param pv_g ...
!> \param pv_nonbond ...
!> \param mm_section ...
!> \author [email protected] (2010-03-01)
!> \note
!> Method: The dipoles are found by minimizing the sum of the electrostatic
!> and the induction energy directly using a conjugate gradient method. This
!> routine assumes that the energy is a quadratic function of the dipoles.
!> Finding the minimum is then done by solving a linear system. This will
!> not work for polarizable force fields that include hyperpolarizability.
!>
!> The implementation of the conjugate gradient solver for linear systems
!> is described in chapter 2.7 Sparse Linear Systems, under the section
!> "Conjugate Gradient Method for a Sparse System". Although the inducible
!> dipoles are the solution of a dense linear system, the same algorithm is
!> still recommended for this situation. One does not have access to the
!> entire hardness kernel to compute the solution with conventional linear
!> algebra routines, but one only has a function that computes the dot
!> product of the hardness kernel and a vector. (This is the routine that
!> computes the electrostatic field at the atoms for a given vector of
!> inducible dipoles.) Given such function, the conjugate gradient method
!> is an efficient way to compute the solution of a linear system, and it
!> scales well towards many degrees of freedom in terms of efficiency and
!> memory usage.
! **************************************************************************************************
SUBROUTINE fist_pol_evaluate_cg(atomic_kind_set, multipoles, ewald_env, ewald_pw, &
fist_nonbond_env, cell, particle_set, local_particles, thermo, vg_coulomb, &
pot_nonbond, f_nonbond, fg_coulomb, use_virial, pv_g, pv_nonbond, mm_section)
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(multipole_type), POINTER :: multipoles
TYPE(ewald_environment_type), POINTER :: ewald_env
TYPE(ewald_pw_type), POINTER :: ewald_pw
TYPE(fist_nonbond_env_type), POINTER :: fist_nonbond_env
TYPE(cell_type), POINTER :: cell
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(distribution_1d_type), POINTER :: local_particles
TYPE(fist_energy_type), POINTER :: thermo
REAL(KIND=dp) :: vg_coulomb, pot_nonbond
REAL(KIND=dp), DIMENSION(:, :) :: f_nonbond, fg_coulomb
LOGICAL, INTENT(IN) :: use_virial
REAL(KIND=dp), DIMENSION(3, 3) :: pv_g, pv_nonbond
TYPE(section_vals_type), POINTER :: mm_section
CHARACTER(len=*), PARAMETER :: routineN = 'fist_pol_evaluate_cg'
INTEGER :: ewald_type, handle, i, iatom, ii, ikind, &
iter, iw, iw2, max_ipol_iter, &
natom_of_kind, natoms, nkind, ntot
INTEGER, DIMENSION(:), POINTER :: atom_list
LOGICAL :: iwarn
REAL(KIND=dp) :: alpha, apol, beta, denom, eps_pol, &
pot_nonbond_local, rmsd
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: conjugate, conjugate_applied, efield1, &
efield1_ext, residual, tmp_dipoles
TYPE(atomic_kind_type), POINTER :: atomic_kind
TYPE(cp_logger_type), POINTER :: logger
CALL timeset(routineN, handle)
NULLIFY (logger, atomic_kind)
logger => cp_get_default_logger()
iw = cp_print_key_unit_nr(logger, mm_section, "PRINT%ITER_INFO", &
extension=".mmLog")
iw2 = cp_print_key_unit_nr(logger, mm_section, "PRINT%EWALD_INFO", &
extension=".mmLog")
CALL ewald_env_get(ewald_env, max_ipol_iter=max_ipol_iter, eps_pol=eps_pol, &
ewald_type=ewald_type)
! allocate work arrays
natoms = SIZE(particle_set)
ALLOCATE (efield1(3, natoms))
ALLOCATE (tmp_dipoles(3, natoms))
ALLOCATE (residual(3, natoms))
ALLOCATE (conjugate(3, natoms))
ALLOCATE (conjugate_applied(3, natoms))
ALLOCATE (efield1_ext(3, natoms))
! Compute the 'external' electrostatic field (all inducible dipoles
! equal to zero). this is required for the conjugate gradient solver.
! We assume that all dipoles are inducible dipoles.
tmp_dipoles(:, :) = multipoles%dipoles ! backup of dipoles
multipoles%dipoles = 0.0_dp
CALL eval_pol_ewald(ewald_type, ewald_env, ewald_pw, fist_nonbond_env, cell, &
particle_set, local_particles, vg_coulomb, pot_nonbond_local, thermo, &
multipoles, do_correction_bonded=.TRUE., do_forces=.FALSE., &
do_stress=.FALSE., do_efield=.TRUE., iw2=iw2, do_debug=.FALSE., &
atomic_kind_set=atomic_kind_set, mm_section=mm_section, &
efield1=efield1_ext)
multipoles%dipoles = tmp_dipoles ! restore backup
! Compute the electric field with the initial guess of the dipoles.
CALL eval_pol_ewald(ewald_type, ewald_env, ewald_pw, fist_nonbond_env, cell, &
particle_set, local_particles, vg_coulomb, pot_nonbond_local, thermo, &
multipoles, do_correction_bonded=.TRUE., do_forces=.FALSE., &
do_stress=.FALSE., do_efield=.TRUE., iw2=iw2, do_debug=.FALSE., &
atomic_kind_set=atomic_kind_set, mm_section=mm_section, &
efield1=efield1)
! Compute the first residual explicitly.
nkind = SIZE(atomic_kind_set)
ntot = 0
residual = 0.0_dp
DO ikind = 1, nkind
atomic_kind => atomic_kind_set(ikind)
CALL get_atomic_kind(atomic_kind, apol=apol, natom=natom_of_kind, atom_list=atom_list)
! ignore atoms with polarizability zero
IF (apol == 0) CYCLE
ntot = ntot + natom_of_kind
DO iatom = 1, natom_of_kind
ii = atom_list(iatom)
DO i = 1, 3
! residual = b - A x
residual(i, ii) = efield1(i, ii) - multipoles%dipoles(i, ii)/apol
END DO
END DO
END DO
! The first conjugate residual is equal to the residual.
conjugate(:, :) = residual
IF (iw > 0) THEN
WRITE (iw, FMT="(/,T2,A,T63,A)") "POL_SCF| Method", "conjugate gradient"
WRITE (iw, FMT="(T2,A,T26,A,T49,A)") "POL_SCF| Iteration", &
"Convergence", "Electrostatic & Induction Energy"
END IF
pol_scf: DO iter = 1, max_ipol_iter
IF (debug_this_module) THEN
! In principle the residual must not be computed explicitly any more. It
! is obtained in an indirect way below. When the DEBUG flag is set, the
! explicit residual is computed and compared with the implicitly derived
! residual as a double check.
CALL eval_pol_ewald(ewald_type, ewald_env, ewald_pw, fist_nonbond_env, cell, &
particle_set, local_particles, vg_coulomb, pot_nonbond_local, thermo, &
multipoles, do_correction_bonded=.TRUE., do_forces=.FALSE., &
do_stress=.FALSE., do_efield=.TRUE., iw2=iw2, do_debug=.FALSE., &
atomic_kind_set=atomic_kind_set, mm_section=mm_section, &
efield1=efield1)
! inapropriate use of denom to check the error on the residual
denom = 0.0_dp
END IF
rmsd = 0.0_dp
! Compute the rmsd of the residual.
DO ikind = 1, nkind
atomic_kind => atomic_kind_set(ikind)
CALL get_atomic_kind(atomic_kind, apol=apol, natom=natom_of_kind, atom_list=atom_list)
! ignore atoms with polarizability zero
IF (apol == 0) CYCLE
DO iatom = 1, natom_of_kind
ii = atom_list(iatom)
DO i = 1, 3
! residual = b - A x
rmsd = rmsd + residual(i, ii)**2
IF (debug_this_module) THEN
denom = denom + (residual(i, ii) - (efield1(i, ii) - &
multipoles%dipoles(i, ii)/apol))**2
END IF
END DO
END DO
END DO
rmsd = SQRT(rmsd/ntot)
IF (iw > 0) THEN
WRITE (iw, FMT="(T2,A,T11,I9,T22,E15.6,T67,A)") "POL_SCF|", iter, rmsd, "(not computed)"
IF (debug_this_module) THEN
denom = SQRT(denom/ntot)
WRITE (iw, FMT="(T2,A,T66,E15.6)") "POL_SCF| Error on implicit residual", denom
END IF
END IF
! Apply the hardness kernel to the conjugate residual.
! We assume that all non-zero dipoles are inducible dipoles.
tmp_dipoles(:, :) = multipoles%dipoles ! backup of dipoles
multipoles%dipoles = conjugate
CALL eval_pol_ewald(ewald_type, ewald_env, ewald_pw, fist_nonbond_env, cell, &
particle_set, local_particles, vg_coulomb, pot_nonbond_local, thermo, &
multipoles, do_correction_bonded=.TRUE., do_forces=.FALSE., &
do_stress=.FALSE., do_efield=.TRUE., iw2=iw2, do_debug=.FALSE., &
atomic_kind_set=atomic_kind_set, mm_section=mm_section, &
efield1=conjugate_applied)
multipoles%dipoles = tmp_dipoles ! restore backup
conjugate_applied(:, :) = efield1_ext - conjugate_applied
! Finish conjugate_applied and compute alpha from the conjugate gradient algorithm.
alpha = 0.0_dp
denom = 0.0_dp
DO ikind = 1, nkind
atomic_kind => atomic_kind_set(ikind)
CALL get_atomic_kind(atomic_kind, apol=apol, natom=natom_of_kind, atom_list=atom_list)
! ignore atoms with polarizability zero
IF (apol == 0) CYCLE
DO iatom = 1, natom_of_kind
ii = atom_list(iatom)
DO i = 1, 3
conjugate_applied(i, ii) = conjugate_applied(i, ii) + conjugate(i, ii)/apol
END DO
alpha = alpha + DOT_PRODUCT(residual(:, ii), residual(:, ii))
denom = denom + DOT_PRODUCT(conjugate(:, ii), conjugate_applied(:, ii))
END DO
END DO
alpha = alpha/denom
! Compute the new residual and beta from the conjugate gradient method.
beta = 0.0_dp
denom = 0.0_dp
DO ikind = 1, nkind
atomic_kind => atomic_kind_set(ikind)
CALL get_atomic_kind(atomic_kind, apol=apol, natom=natom_of_kind, atom_list=atom_list)
IF (apol == 0) CYCLE
DO iatom = 1, natom_of_kind
ii = atom_list(iatom)
denom = denom + DOT_PRODUCT(residual(:, ii), residual(:, ii))
DO i = 1, 3
residual(i, ii) = residual(i, ii) - alpha*conjugate_applied(i, ii)
END DO
beta = beta + DOT_PRODUCT(residual(:, ii), residual(:, ii))
END DO
END DO
beta = beta/denom
! Compute the new dipoles, the new conjugate residual, and the induction
! energy.
thermo%e_induction = 0.0_dp
DO ikind = 1, nkind
atomic_kind => atomic_kind_set(ikind)
CALL get_atomic_kind(atomic_kind, apol=apol, natom=natom_of_kind, atom_list=atom_list)
! ignore atoms with polarizability zero
IF (apol == 0) CYCLE
DO iatom = 1, natom_of_kind
ii = atom_list(iatom)
DO i = 1, 3
multipoles%dipoles(i, ii) = multipoles%dipoles(i, ii) + alpha*conjugate(i, ii)
conjugate(i, ii) = residual(i, ii) + beta*conjugate(i, ii)
thermo%e_induction = thermo%e_induction + multipoles%dipoles(i, ii)**2/apol/2.0_dp
END DO
END DO
END DO
! Quit if rmsd is low enough
IF (rmsd <= eps_pol) THEN
IF (iw > 0) WRITE (iw, FMT="(T2,A)") "POL_SCF| Self-consistent polarization converged"
EXIT pol_scf
END IF
! Print warning when not converged
iwarn = ((rmsd > eps_pol) .AND. (iter >= max_ipol_iter))
IF (iwarn) THEN
IF (iw > 0) THEN
WRITE (iw, FMT="(T2,A,I0,A,ES9.3)") &
"POL_SCF| Self-consistent polarization not converged in ", max_ipol_iter, &
" steps to ", eps_pol
END IF
CPWARN("Self-consistent Polarization not converged!")
END IF
END DO pol_scf
IF (debug_this_module) THEN
! Now evaluate after convergence to obtain forces and converged energies
CALL eval_pol_ewald(ewald_type, ewald_env, ewald_pw, fist_nonbond_env, cell, &
particle_set, local_particles, vg_coulomb, pot_nonbond_local, thermo, &
multipoles, do_correction_bonded=.TRUE., do_forces=.TRUE., &
do_stress=use_virial, do_efield=.TRUE., iw2=iw2, do_debug=.FALSE., &
atomic_kind_set=atomic_kind_set, mm_section=mm_section, &
forces_local=fg_coulomb, forces_glob=f_nonbond, &
pv_local=pv_g, pv_glob=pv_nonbond, efield1=efield1)
! Do a final check on the convergence: compute the residual explicitely
rmsd = 0.0_dp
DO ikind = 1, nkind
atomic_kind => atomic_kind_set(ikind)
CALL get_atomic_kind(atomic_kind, apol=apol, natom=natom_of_kind, atom_list=atom_list)
! ignore atoms with polarizability zero
IF (apol == 0) CYCLE
DO iatom = 1, natom_of_kind
ii = atom_list(iatom)
DO i = 1, 3
! residual = b - A x
rmsd = rmsd + (efield1(i, ii) - multipoles%dipoles(i, ii)/apol)**2
END DO
END DO
END DO
rmsd = SQRT(rmsd/ntot)
IF (iw > 0) WRITE (iw, FMT="(T2,A,T66,E15.6)") "POL_SCF| Final RMSD of residual", rmsd
! Stop program when congergence is not reached after all
IF (rmsd > eps_pol) THEN
CPABORT("Error in the conjugate gradient method for self-consistent polarization!")
END IF
ELSE
! Now evaluate after convergence to obtain forces and converged energies
CALL eval_pol_ewald(ewald_type, ewald_env, ewald_pw, fist_nonbond_env, cell, &
particle_set, local_particles, vg_coulomb, pot_nonbond_local, thermo, &
multipoles, do_correction_bonded=.TRUE., do_forces=.TRUE., &
do_stress=use_virial, do_efield=.FALSE., iw2=iw2, do_debug=.FALSE., &
atomic_kind_set=atomic_kind_set, mm_section=mm_section, &
forces_local=fg_coulomb, forces_glob=f_nonbond, &
pv_local=pv_g, pv_glob=pv_nonbond)
END IF
pot_nonbond = pot_nonbond + pot_nonbond_local
CALL logger%para_env%sum(pot_nonbond_local)
IF (iw > 0) WRITE (iw, FMT="(T2,A,T61,F20.10)") "POL_SCF| Final", &
vg_coulomb + pot_nonbond_local + thermo%e_induction
! Deallocate working arrays
DEALLOCATE (efield1)
DEALLOCATE (tmp_dipoles)
DEALLOCATE (residual)
DEALLOCATE (conjugate)
DEALLOCATE (conjugate_applied)
DEALLOCATE (efield1_ext)
CALL cp_print_key_finished_output(iw2, logger, mm_section, &
"PRINT%EWALD_INFO")
CALL cp_print_key_finished_output(iw, logger, mm_section, &
"PRINT%ITER_INFO")
CALL timestop(handle)
END SUBROUTINE fist_pol_evaluate_cg
! **************************************************************************************************
!> \brief Main driver for evaluating electrostatic in polarible forcefields
!> All the dependence on the Ewald method should go here!
!> \param ewald_type ...
!> \param ewald_env ...
!> \param ewald_pw ...
!> \param fist_nonbond_env ...
!> \param cell ...
!> \param particle_set ...
!> \param local_particles ...
!> \param vg_coulomb ...
!> \param pot_nonbond ...
!> \param thermo ...
!> \param multipoles ...
!> \param do_correction_bonded ...
!> \param do_forces ...
!> \param do_stress ...
!> \param do_efield ...
!> \param iw2 ...
!> \param do_debug ...
!> \param atomic_kind_set ...
!> \param mm_section ...
!> \param efield0 ...
!> \param efield1 ...
!> \param efield2 ...
!> \param forces_local ...
!> \param forces_glob ...
!> \param pv_local ...
!> \param pv_glob ...
!> \author Teodoro Laino [tlaino] 05.2009
! **************************************************************************************************
SUBROUTINE eval_pol_ewald(ewald_type, ewald_env, ewald_pw, fist_nonbond_env, &
cell, particle_set, local_particles, vg_coulomb, pot_nonbond, thermo, &
multipoles, do_correction_bonded, do_forces, do_stress, do_efield, iw2, &
do_debug, atomic_kind_set, mm_section, efield0, efield1, efield2, forces_local, &
forces_glob, pv_local, pv_glob)
INTEGER, INTENT(IN) :: ewald_type
TYPE(ewald_environment_type), POINTER :: ewald_env
TYPE(ewald_pw_type), POINTER :: ewald_pw
TYPE(fist_nonbond_env_type), POINTER :: fist_nonbond_env
TYPE(cell_type), POINTER :: cell
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(distribution_1d_type), POINTER :: local_particles
REAL(KIND=dp), INTENT(OUT) :: vg_coulomb, pot_nonbond
TYPE(fist_energy_type), POINTER :: thermo
TYPE(multipole_type), POINTER :: multipoles
LOGICAL, INTENT(IN) :: do_correction_bonded, do_forces, &
do_stress, do_efield
INTEGER, INTENT(IN) :: iw2
LOGICAL, INTENT(IN) :: do_debug
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(section_vals_type), POINTER :: mm_section
REAL(KIND=dp), DIMENSION(:), OPTIONAL :: efield0
REAL(KIND=dp), DIMENSION(:, :), OPTIONAL :: efield1, efield2
REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT), &
OPTIONAL :: forces_local, forces_glob, pv_local, &
pv_glob
CHARACTER(len=*), PARAMETER :: routineN = 'eval_pol_ewald'
INTEGER :: handle
CALL timeset(routineN, handle)
pot_nonbond = 0.0_dp ! Initialization..
vg_coulomb = 0.0_dp ! Initialization..
SELECT CASE (ewald_type)
CASE (do_ewald_ewald)
CALL ewald_multipole_evaluate(ewald_env, ewald_pw, fist_nonbond_env, cell, &
particle_set, local_particles, vg_coulomb, pot_nonbond, thermo%e_neut, &
thermo%e_self, multipoles%task, do_correction_bonded=do_correction_bonded, &
do_forces=do_forces, do_stress=do_stress, do_efield=do_efield, &
radii=multipoles%radii, charges=multipoles%charges, &
dipoles=multipoles%dipoles, quadrupoles=multipoles%quadrupoles, &
forces_local=forces_local, forces_glob=forces_glob, pv_local=pv_local, &
pv_glob=pv_glob, iw=iw2, do_debug=do_debug, atomic_kind_set=atomic_kind_set, &
mm_section=mm_section, efield0=efield0, efield1=efield1, efield2=efield2)
CASE (do_ewald_pme)
CPABORT("Multipole Ewald not yet implemented within a PME scheme!")
CASE (do_ewald_spme)
CPABORT("Multipole Ewald not yet implemented within a SPME scheme!")
END SELECT
CALL timestop(handle)
END SUBROUTINE eval_pol_ewald
END MODULE fist_pol_scf