-
Notifications
You must be signed in to change notification settings - Fork 1
/
fist_force.F
906 lines (844 loc) · 45.1 KB
/
fist_force.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
!--------------------------------------------------------------------------------------------------!
! 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
!> Torsions added(DG)05-Dec-2000
!> Variable names changed(DG)05-Dec-2000
!> CJM SEPT-12-2002: int_env is now passed
!> CJM NOV-30-2003: only uses fist_env
!> \author CJM & JGH
! **************************************************************************************************
MODULE fist_force
USE atomic_kind_types, ONLY: atomic_kind_type,&
get_atomic_kind
USE atprop_types, ONLY: atprop_type
USE cell_methods, ONLY: init_cell
USE cell_types, ONLY: cell_type
USE cp_log_handling, ONLY: cp_get_default_logger,&
cp_logger_type
USE cp_output_handling, ONLY: cp_p_file,&
cp_print_key_finished_output,&
cp_print_key_should_output,&
cp_print_key_unit_nr
USE cp_subsys_types, ONLY: cp_subsys_get,&
cp_subsys_type
USE cp_units, ONLY: cp_unit_from_cp2k
USE distribution_1d_types, ONLY: distribution_1d_type
USE ewald_environment_types, ONLY: ewald_env_get,&
ewald_environment_type
USE ewald_pw_methods, ONLY: ewald_pw_grid_update
USE ewald_pw_types, ONLY: ewald_pw_type
USE ewalds, ONLY: ewald_evaluate,&
ewald_print,&
ewald_self,&
ewald_self_atom
USE ewalds_multipole, ONLY: ewald_multipole_evaluate
USE exclusion_types, ONLY: exclusion_type
USE fist_efield_methods, ONLY: fist_dipole,&
fist_efield_energy_force
USE fist_efield_types, ONLY: fist_efield_type
USE fist_energy_types, ONLY: fist_energy_type
USE fist_environment_types, ONLY: fist_env_get,&
fist_environment_type
USE fist_intra_force, ONLY: force_intra_control
USE fist_neighbor_list_control, ONLY: list_control
USE fist_nonbond_env_types, ONLY: fist_nonbond_env_type
USE fist_nonbond_force, ONLY: bonded_correct_gaussian,&
force_nonbond
USE fist_pol_scf, ONLY: fist_pol_evaluate
USE input_constants, ONLY: do_fist_pol_none
USE input_section_types, ONLY: section_vals_get_subs_vals,&
section_vals_type
USE kinds, ONLY: dp
USE manybody_eam, ONLY: density_nonbond
USE manybody_potential, ONLY: energy_manybody,&
force_nonbond_manybody
USE message_passing, ONLY: mp_para_env_type
USE molecule_kind_types, ONLY: molecule_kind_type
USE molecule_types, ONLY: molecule_type
USE multipole_types, ONLY: multipole_type
USE particle_types, ONLY: particle_type
USE pme, ONLY: pme_evaluate
USE pw_poisson_types, ONLY: do_ewald_ewald,&
do_ewald_none,&
do_ewald_pme,&
do_ewald_spme
USE shell_potential_types, ONLY: shell_kind_type
USE spme, ONLY: spme_evaluate
USE virial_types, ONLY: virial_type
#include "./base/base_uses.f90"
IMPLICIT NONE
PRIVATE
PUBLIC :: fist_calc_energy_force
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'fist_force'
TYPE debug_variables_type
REAL(KIND=dp) :: pot_manybody = 0.0_dp, pot_bend = 0.0_dp, pot_torsion = 0.0_dp
REAL(KIND=dp) :: pot_nonbond = 0.0_dp, pot_g = 0.0_dp, pot_bond = 0.0_dp
REAL(KIND=dp) :: pot_imptors = 0.0_dp, pot_urey_bradley = 0.0_dp
REAL(KIND=dp) :: pot_opbend = 0.0_dp
REAL(KIND=dp), DIMENSION(:, :), POINTER :: f_nonbond => NULL(), f_g => NULL(), f_bond => NULL(), f_bend => NULL(), &
f_torsion => NULL(), f_imptors => NULL(), f_ub => NULL(), &
f_opbend => NULL()
REAL(KIND=dp), DIMENSION(3, 3) :: pv_nonbond = 0.0_dp, pv_g = 0.0_dp, pv_bond = 0.0_dp, pv_ub = 0.0_dp, &
pv_bend = 0.0_dp, pv_torsion = 0.0_dp, pv_imptors = 0.0_dp, &
pv_opbend = 0.0_dp
END TYPE debug_variables_type
CONTAINS
! **************************************************************************************************
!> \brief Calculates the total potential energy, total force, and the
!> total pressure tensor from the potentials
!> \param fist_env ...
!> \param debug ...
!> \par History
!> Harald Forbert(Dec-2000): Changes for multiple linked lists
!> cjm, 20-Feb-2001: box_ref used to initialize ewald. Now
!> have consistent restarts with npt and ewald
!> JGH(15-Mar-2001): box_change replaces ensemble parameter
!> Call ewald_setup if box has changed
!> Consistent setup for PME and SPME
!> cjm, 28-Feb-2006: box_change is gone
!> \author CJM & JGH
! **************************************************************************************************
SUBROUTINE fist_calc_energy_force(fist_env, debug)
TYPE(fist_environment_type), POINTER :: fist_env
TYPE(debug_variables_type), INTENT(INOUT), &
OPTIONAL :: debug
CHARACTER(len=*), PARAMETER :: routineN = 'fist_calc_energy_force'
INTEGER :: do_ipol, ewald_type, fg_coulomb_size, handle, i, ii, ikind, iparticle_kind, &
iparticle_local, iw, iw2, j, natoms, nlocal_particles, node, nparticle_kind, &
nparticle_local, nshell, shell_index
LOGICAL :: do_multipoles, shell_model_ad, &
shell_present, use_virial
REAL(KIND=dp) :: ef_ener, fc, fs, mass, pot_bend, pot_bond, pot_imptors, pot_manybody, &
pot_nonbond, pot_opbend, pot_shell, pot_torsion, pot_urey_bradley, vg_coulomb, xdum1, &
xdum2, xdum3
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: ef_f, f_nonbond, f_total, fcore_nonbond, &
fcore_shell_bonded, fcore_total, fg_coulomb, fgcore_coulomb, fgshell_coulomb, &
fshell_nonbond, fshell_total
REAL(KIND=dp), DIMENSION(3, 3) :: ef_pv, ident, pv_bc, pv_bend, pv_bond, &
pv_g, pv_imptors, pv_nonbond, &
pv_opbend, pv_torsion, pv_urey_bradley
REAL(KIND=dp), DIMENSION(:), POINTER :: e_coulomb
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(atomic_kind_type), POINTER :: atomic_kind
TYPE(atprop_type), POINTER :: atprop_env
TYPE(cell_type), POINTER :: cell
TYPE(cp_logger_type), POINTER :: logger
TYPE(cp_subsys_type), POINTER :: subsys
TYPE(distribution_1d_type), POINTER :: local_molecules, local_particles
TYPE(ewald_environment_type), POINTER :: ewald_env
TYPE(ewald_pw_type), POINTER :: ewald_pw
TYPE(exclusion_type), DIMENSION(:), POINTER :: exclusions
TYPE(fist_efield_type), POINTER :: efield
TYPE(fist_energy_type), POINTER :: thermo
TYPE(fist_nonbond_env_type), POINTER :: fist_nonbond_env
TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
TYPE(mp_para_env_type), POINTER :: para_env
TYPE(multipole_type), POINTER :: multipoles
TYPE(particle_type), DIMENSION(:), POINTER :: core_particle_set, particle_set, &
shell_particle_set
TYPE(section_vals_type), POINTER :: force_env_section, mm_section, &
print_section
TYPE(shell_kind_type), POINTER :: shell
TYPE(virial_type), POINTER :: virial
CALL timeset(routineN, handle)
NULLIFY (logger)
NULLIFY (subsys, virial, atprop_env, para_env, force_env_section)
logger => cp_get_default_logger()
CALL fist_env_get(fist_env, &
subsys=subsys, &
para_env=para_env, &
input=force_env_section)
CALL cp_subsys_get(subsys, &
virial=virial, &
atprop=atprop_env)
use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
NULLIFY (atomic_kind, atomic_kind_set, cell, ewald_pw, ewald_env, &
fist_nonbond_env, mm_section, local_molecules, local_particles, &
molecule_kind_set, molecule_set, particle_set, print_section, &
shell, shell_particle_set, core_particle_set, thermo, multipoles, &
e_coulomb)
mm_section => section_vals_get_subs_vals(force_env_section, "MM")
iw = cp_print_key_unit_nr(logger, mm_section, "PRINT%DERIVATIVES", &
extension=".mmLog")
iw2 = cp_print_key_unit_nr(logger, mm_section, "PRINT%EWALD_INFO", &
extension=".mmLog")
CALL fist_env_get(fist_env, ewald_pw=ewald_pw, ewald_env=ewald_env, &
local_particles=local_particles, particle_set=particle_set, &
atomic_kind_set=atomic_kind_set, molecule_set=molecule_set, &
local_molecules=local_molecules, thermo=thermo, &
molecule_kind_set=molecule_kind_set, fist_nonbond_env=fist_nonbond_env, &
cell=cell, shell_model=shell_present, shell_model_ad=shell_model_ad, &
multipoles=multipoles, exclusions=exclusions, efield=efield)
CALL ewald_env_get(ewald_env, ewald_type=ewald_type, do_multipoles=do_multipoles, &
do_ipol=do_ipol)
! Initialize ewald grids
CALL init_cell(cell)
CALL ewald_pw_grid_update(ewald_pw, ewald_env, cell%hmat)
natoms = SIZE(particle_set)
nlocal_particles = 0
nparticle_kind = SIZE(atomic_kind_set)
DO ikind = 1, nparticle_kind
nlocal_particles = nlocal_particles + local_particles%n_el(ikind)
END DO
ALLOCATE (f_nonbond(3, natoms))
f_nonbond = 0.0_dp
nshell = 0
IF (shell_present) THEN
CALL fist_env_get(fist_env, shell_particle_set=shell_particle_set, &
core_particle_set=core_particle_set)
CPASSERT(ASSOCIATED(shell_particle_set))
CPASSERT(ASSOCIATED(core_particle_set))
nshell = SIZE(shell_particle_set)
ALLOCATE (fshell_nonbond(3, nshell))
fshell_nonbond = 0.0_dp
ALLOCATE (fcore_nonbond(3, nshell))
fcore_nonbond = 0.0_dp
ELSE
NULLIFY (shell_particle_set, core_particle_set)
END IF
IF (fist_nonbond_env%do_nonbonded) THEN
! First check with list_control to update neighbor lists
IF (ASSOCIATED(exclusions)) THEN
CALL list_control(atomic_kind_set, particle_set, local_particles, &
cell, fist_nonbond_env, para_env, mm_section, shell_particle_set, &
core_particle_set, exclusions=exclusions)
ELSE
CALL list_control(atomic_kind_set, particle_set, local_particles, &
cell, fist_nonbond_env, para_env, mm_section, shell_particle_set, &
core_particle_set)
END IF
END IF
! Initialize force, energy and pressure tensor arrays
DO i = 1, natoms
particle_set(i)%f(1) = 0.0_dp
particle_set(i)%f(2) = 0.0_dp
particle_set(i)%f(3) = 0.0_dp
END DO
IF (nshell > 0) THEN
DO i = 1, nshell
shell_particle_set(i)%f(1) = 0.0_dp
shell_particle_set(i)%f(2) = 0.0_dp
shell_particle_set(i)%f(3) = 0.0_dp
core_particle_set(i)%f(1) = 0.0_dp
core_particle_set(i)%f(2) = 0.0_dp
core_particle_set(i)%f(3) = 0.0_dp
END DO
END IF
IF (use_virial) THEN
pv_bc = 0.0_dp
pv_bond = 0.0_dp
pv_bend = 0.0_dp
pv_torsion = 0.0_dp
pv_imptors = 0.0_dp
pv_opbend = 0.0_dp
pv_urey_bradley = 0.0_dp
pv_nonbond = 0.0_dp
pv_g = 0.0_dp
virial%pv_virial = 0.0_dp
END IF
pot_nonbond = 0.0_dp
pot_manybody = 0.0_dp
pot_bond = 0.0_dp
pot_bend = 0.0_dp
pot_torsion = 0.0_dp
pot_opbend = 0.0_dp
pot_imptors = 0.0_dp
pot_urey_bradley = 0.0_dp
pot_shell = 0.0_dp
vg_coulomb = 0.0_dp
thermo%pot = 0.0_dp
thermo%harm_shell = 0.0_dp
! Get real-space non-bonded forces:
IF (iw > 0) THEN
WRITE (iw, '(A)') " FIST:: FORCES IN INPUT..."
WRITE (iw, '(3f15.9)') ((particle_set(i)%f(j), j=1, 3), i=1, SIZE(particle_set))
END IF
IF (fist_nonbond_env%do_nonbonded) THEN
! Compute density for EAM
CALL density_nonbond(fist_nonbond_env, particle_set, cell, para_env)
! Compute embedding function and manybody energy
CALL energy_manybody(fist_nonbond_env, atomic_kind_set, local_particles, particle_set, &
cell, pot_manybody, para_env, mm_section, use_virial)
! Nonbond contribution + Manybody Forces
IF (shell_present) THEN
CALL force_nonbond(fist_nonbond_env, ewald_env, particle_set, cell, &
pot_nonbond, f_nonbond, pv_nonbond, &
fshell_nonbond=fshell_nonbond, fcore_nonbond=fcore_nonbond, &
atprop_env=atprop_env, &
atomic_kind_set=atomic_kind_set, use_virial=use_virial)
ELSE
CALL force_nonbond(fist_nonbond_env, ewald_env, particle_set, cell, &
pot_nonbond, f_nonbond, pv_nonbond, atprop_env=atprop_env, &
atomic_kind_set=atomic_kind_set, use_virial=use_virial)
CALL force_nonbond_manybody(fist_nonbond_env, particle_set, cell, f_nonbond, pv_nonbond, &
use_virial=use_virial)
END IF
END IF
IF (iw > 0) THEN
WRITE (iw, '(A)') " FIST:: NONBOND + R-SPACE ELECTROSTATIC FORCES ..."
WRITE (iw, '(3f15.9)') f_nonbond
IF (shell_present .AND. shell_model_ad) THEN
WRITE (iw, '(A)') " FIST:: SHELL NONBOND + R-SPACE ELECTROSTATIC FORCES ..."
WRITE (iw, '(3f15.9)') fshell_nonbond
WRITE (iw, '(A)') " FIST:: CORE NONBOND + R-SPACE ELECTROSTATIC FORCES ..."
WRITE (iw, '(3f15.9)') fcore_nonbond
END IF
END IF
! Get g-space non-bonded forces:
IF (ewald_type /= do_ewald_none) THEN
! Determine size of the forces array
SELECT CASE (ewald_type)
CASE (do_ewald_ewald)
fg_coulomb_size = nlocal_particles
CASE DEFAULT
fg_coulomb_size = natoms
END SELECT
! Allocate and zeroing arrays
ALLOCATE (fg_coulomb(3, fg_coulomb_size))
fg_coulomb = 0.0_dp
IF (shell_present) THEN
ALLOCATE (fgshell_coulomb(3, nshell))
ALLOCATE (fgcore_coulomb(3, nshell))
fgshell_coulomb = 0.0_dp
fgcore_coulomb = 0.0_dp
END IF
IF (shell_present .AND. do_multipoles) THEN
CPABORT("Multipoles and Core-Shell model not implemented.")
END IF
! If not multipole: Compute self-interaction and neutralizing background
! for multipoles is handled separately..
IF (.NOT. do_multipoles) THEN
CALL ewald_self(ewald_env, cell, atomic_kind_set, local_particles, &
thermo%e_self, thermo%e_neut, fist_nonbond_env%charges)
IF (atprop_env%energy) THEN
CALL ewald_self_atom(ewald_env, atomic_kind_set, local_particles, &
atprop_env%atener, fist_nonbond_env%charges)
atprop_env%atener = atprop_env%atener + thermo%e_neut/SIZE(atprop_env%atener)
END IF
END IF
! Polarizable force-field
IF (do_ipol /= do_fist_pol_none) THEN
! Check if an array of chagres was provided and in case abort due to lack of implementation
IF (ASSOCIATED(fist_nonbond_env%charges)) THEN
CPABORT("Polarizable force field and array charges not implemented!")
END IF
! Converge the dipoles self-consistently
CALL 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)
ELSE
! Non-Polarizable force-field
SELECT CASE (ewald_type)
CASE (do_ewald_ewald)
! Parallelisation over atoms --> allocate local atoms
IF (shell_present) THEN
! Check if an array of chagres was provided and in case abort due to lack of implementation
IF (ASSOCIATED(fist_nonbond_env%charges)) THEN
CPABORT("Core-Shell and array charges not implemented!")
END IF
IF (do_multipoles) THEN
CPABORT("Multipole Ewald and CORE-SHELL not yet implemented within Ewald sum!")
ELSE
CPABORT("Core-Shell model not yet implemented within Ewald sums.")
END IF
ELSE
IF (do_multipoles) THEN
! Check if an array of chagres was provided and in case abort due to lack of implementation
IF (ASSOCIATED(fist_nonbond_env%charges)) THEN
CPABORT("Multipole Ewald and array charges not implemented!")
END IF
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=.TRUE., do_forces=.TRUE., &
do_stress=use_virial, do_efield=.FALSE., radii=multipoles%radii, &
charges=multipoles%charges, dipoles=multipoles%dipoles, &
quadrupoles=multipoles%quadrupoles, forces_local=fg_coulomb, &
forces_glob=f_nonbond, pv_local=pv_g, pv_glob=pv_nonbond, iw=iw2, &
do_debug=.TRUE., atomic_kind_set=atomic_kind_set, mm_section=mm_section)
ELSE
IF (atprop_env%energy) THEN
ALLOCATE (e_coulomb(fg_coulomb_size))
END IF
CALL ewald_evaluate(ewald_env, ewald_pw, cell, atomic_kind_set, particle_set, &
local_particles, fg_coulomb, vg_coulomb, pv_g, use_virial=use_virial, &
charges=fist_nonbond_env%charges, e_coulomb=e_coulomb)
END IF
END IF
CASE (do_ewald_pme)
! Parallelisation over grids --> allocate all atoms
IF (shell_present) THEN
! Check if an array of chagres was provided and in case abort due to lack of implementation
IF (ASSOCIATED(fist_nonbond_env%charges)) THEN
CPABORT("Core-Shell and array charges not implemented!")
END IF
IF (do_multipoles) THEN
CPABORT("Multipole Ewald and CORE-SHELL not yet implemented within a PME scheme!")
ELSE
CALL pme_evaluate(ewald_env, ewald_pw, cell, particle_set, vg_coulomb, fg_coulomb, &
pv_g, shell_particle_set=shell_particle_set, core_particle_set=core_particle_set, &
fgshell_coulomb=fgshell_coulomb, fgcore_coulomb=fgcore_coulomb, use_virial=use_virial, &
atprop=atprop_env)
CALL para_env%sum(fgshell_coulomb)
CALL para_env%sum(fgcore_coulomb)
END IF
ELSE
IF (do_multipoles) THEN
CPABORT("Multipole Ewald not yet implemented within a PME scheme!")
ELSE
CALL pme_evaluate(ewald_env, ewald_pw, cell, particle_set, vg_coulomb, fg_coulomb, &
pv_g, use_virial=use_virial, charges=fist_nonbond_env%charges, &
atprop=atprop_env)
END IF
END IF
CALL para_env%sum(fg_coulomb)
CASE (do_ewald_spme)
! Parallelisation over grids --> allocate all atoms
IF (shell_present) THEN
! Check if an array of charges was provided and in case abort due to lack of implementation
IF (ASSOCIATED(fist_nonbond_env%charges)) THEN
CPABORT("Core-Shell and array charges not implemented!")
END IF
IF (do_multipoles) THEN
CPABORT("Multipole Ewald and CORE-SHELL not yet implemented within a SPME scheme!")
ELSE
CALL spme_evaluate(ewald_env, ewald_pw, cell, particle_set, fg_coulomb, vg_coulomb, &
pv_g, shell_particle_set=shell_particle_set, core_particle_set=core_particle_set, &
fgshell_coulomb=fgshell_coulomb, fgcore_coulomb=fgcore_coulomb, use_virial=use_virial, &
atprop=atprop_env)
CALL para_env%sum(fgshell_coulomb)
CALL para_env%sum(fgcore_coulomb)
END IF
ELSE
IF (do_multipoles) THEN
CPABORT("Multipole Ewald not yet implemented within a SPME scheme!")
ELSE
CALL spme_evaluate(ewald_env, ewald_pw, cell, particle_set, fg_coulomb, vg_coulomb, &
pv_g, use_virial=use_virial, charges=fist_nonbond_env%charges, &
atprop=atprop_env)
END IF
END IF
CALL para_env%sum(fg_coulomb)
END SELECT
END IF
! Subtract the interaction between screening charges. This is a
! correction in real-space and depends on the neighborlists. Therefore
! it is only executed if fist_nonbond_env%do_nonbonded is set.
IF (fist_nonbond_env%do_nonbonded) THEN
! Correction for core-shell model
IF (shell_present) THEN
CALL bonded_correct_gaussian(fist_nonbond_env, atomic_kind_set, &
local_particles, particle_set, ewald_env, thermo%e_bonded, &
pv_bc, shell_particle_set=shell_particle_set, &
core_particle_set=core_particle_set, atprop_env=atprop_env, cell=cell, &
use_virial=use_virial)
ELSE
IF (.NOT. do_multipoles) THEN
CALL bonded_correct_gaussian(fist_nonbond_env, &
atomic_kind_set, local_particles, particle_set, &
ewald_env, thermo%e_bonded, pv_bc=pv_bc, atprop_env=atprop_env, cell=cell, &
use_virial=use_virial)
END IF
END IF
END IF
IF (.NOT. do_multipoles) THEN
! Multipole code has its own printing routine.
CALL ewald_print(iw2, pot_nonbond, vg_coulomb, thermo%e_self, thermo%e_neut, &
thermo%e_bonded)
END IF
ELSE
IF (use_virial) THEN
pv_g = 0.0_dp
pv_bc = 0.0_dp
END IF
thermo%e_neut = 0.0_dp
END IF
IF (iw > 0) THEN
IF (ALLOCATED(fg_coulomb)) THEN
WRITE (iw, '(A)') " FIST:: NONBONDED ELECTROSTATIC FORCES IN G-SPACE..."
WRITE (iw, '(3f15.9)') ((fg_coulomb(j, i), j=1, 3), i=1, SIZE(fg_coulomb, 2))
END IF
IF (shell_present .AND. shell_model_ad .AND. ALLOCATED(fgshell_coulomb)) THEN
WRITE (iw, '(A)') " FIST:: SHELL NONBONDED ELECTROSTATIC FORCES IN G-SPACE..."
WRITE (iw, '(3f15.9)') ((fgshell_coulomb(j, i), j=1, 3), i=1, SIZE(fgshell_coulomb, 2))
WRITE (iw, '(A)') " FIST:: CORE NONBONDED ELECTROSTATIC FORCES IN G-SPACE..."
WRITE (iw, '(3f15.9)') ((fgcore_coulomb(j, i), j=1, 3), i=1, SIZE(fgcore_coulomb, 2))
END IF
END IF
! calculate the action of an (external) electric field
IF (efield%apply_field) THEN
ALLOCATE (ef_f(3, natoms))
CALL fist_efield_energy_force(ef_ener, ef_f, ef_pv, atomic_kind_set, particle_set, cell, &
efield, use_virial=use_virial, iunit=iw, charges=fist_nonbond_env%charges)
END IF
! Get intramolecular forces
IF (PRESENT(debug)) THEN
CALL force_intra_control(molecule_set, molecule_kind_set, &
local_molecules, particle_set, shell_particle_set, &
core_particle_set, pot_bond, pot_bend, pot_urey_bradley, &
pot_torsion, pot_imptors, pot_opbend, pot_shell, pv_bond, pv_bend, &
pv_urey_bradley, pv_torsion, pv_imptors, pv_opbend, &
debug%f_bond, debug%f_bend, debug%f_torsion, debug%f_ub, &
debug%f_imptors, debug%f_opbend, cell, use_virial, atprop_env)
ELSE
CALL force_intra_control(molecule_set, molecule_kind_set, &
local_molecules, particle_set, shell_particle_set, &
core_particle_set, pot_bond, pot_bend, pot_urey_bradley, &
pot_torsion, pot_imptors, pot_opbend, pot_shell, pv_bond, pv_bend, &
pv_urey_bradley, pv_torsion, pv_imptors, pv_opbend, &
cell=cell, use_virial=use_virial, atprop_env=atprop_env)
END IF
! Perform global sum of the intra-molecular (bonded) forces for the core-shell atoms
IF (shell_present .AND. shell_model_ad) THEN
ALLOCATE (fcore_shell_bonded(3, nshell))
fcore_shell_bonded = 0.0_dp
DO i = 1, natoms
shell_index = particle_set(i)%shell_index
IF (shell_index /= 0) THEN
fcore_shell_bonded(1:3, shell_index) = particle_set(i)%f(1:3)
END IF
END DO
CALL para_env%sum(fcore_shell_bonded)
DO i = 1, natoms
shell_index = particle_set(i)%shell_index
IF (shell_index /= 0) THEN
particle_set(i)%f(1:3) = fcore_shell_bonded(1:3, shell_index)
END IF
END DO
DEALLOCATE (fcore_shell_bonded)
END IF
IF (iw > 0) THEN
xdum1 = cp_unit_from_cp2k(pot_bond, "kcalmol")
xdum2 = cp_unit_from_cp2k(pot_bend, "kcalmol")
xdum3 = cp_unit_from_cp2k(pot_urey_bradley, "kcalmol")
WRITE (iw, '(A)') " FIST energy contributions in kcal/mol:"
WRITE (iw, '(1x,"BOND = ",f13.4,'// &
'2x,"ANGLE = ",f13.4,'// &
'2x,"UBRAD = ",f13.4)') xdum1, xdum2, xdum3
xdum1 = cp_unit_from_cp2k(pot_torsion, "kcalmol")
xdum2 = cp_unit_from_cp2k(pot_imptors, "kcalmol")
xdum3 = cp_unit_from_cp2k(pot_opbend, "kcalmol")
WRITE (iw, '(1x,"TORSION = ",f13.4,'// &
'2x,"IMPTORS = ",f13.4,'// &
'2x,"OPBEND = ",f13.4)') xdum1, xdum2, xdum3
WRITE (iw, '(A)') " FIST:: CORRECTED BONDED ELECTROSTATIC FORCES + INTERNAL FORCES..."
WRITE (iw, '(3f15.9)') ((particle_set(i)%f(j), j=1, 3), i=1, SIZE(particle_set))
IF (shell_present .AND. shell_model_ad) THEN
WRITE (iw, '(A)') " FIST:: SHELL CORRECTED BONDED ELECTROSTATIC FORCES + INTERNAL FORCES..."
WRITE (iw, '(3f15.9)') ((shell_particle_set(i)%f(j), j=1, 3), i=1, SIZE(shell_particle_set))
WRITE (iw, '(A)') " FIST:: CORE CORRECTED BONDED ELECTROSTATIC FORCES + INTERNAL FORCES..."
WRITE (iw, '(3f15.9)') ((core_particle_set(i)%f(j), j=1, 3), i=1, SIZE(core_particle_set))
END IF
END IF
! add up all the potential energies
thermo%pot = pot_nonbond + pot_bond + pot_bend + pot_torsion + pot_opbend + &
pot_imptors + pot_urey_bradley + pot_manybody + pot_shell
CALL para_env%sum(thermo%pot)
IF (shell_present) THEN
thermo%harm_shell = pot_shell
CALL para_env%sum(thermo%harm_shell)
END IF
! add g-space contributions if needed
IF (ewald_type /= do_ewald_none) THEN
! e_self, e_neut, and ebonded are already summed over all processors
! vg_coulomb is not calculated in parallel
thermo%e_gspace = vg_coulomb
thermo%pot = thermo%pot + thermo%e_self + thermo%e_neut
thermo%pot = thermo%pot + vg_coulomb + thermo%e_bonded
! add the induction energy of the dipoles for polarizable models
IF (do_ipol /= do_fist_pol_none) thermo%pot = thermo%pot + thermo%e_induction
END IF
! add up all the forces
!
! nonbonded forces might be calculated for atoms not on this node
! ewald forces are strictly local -> sum only over pnode
! We first sum the forces in f_nonbond, this allows for a more efficient
! global sum in the parallel code and in the end copy them back to part
ALLOCATE (f_total(3, natoms))
f_total = 0.0_dp
DO i = 1, natoms
f_total(1, i) = particle_set(i)%f(1) + f_nonbond(1, i)
f_total(2, i) = particle_set(i)%f(2) + f_nonbond(2, i)
f_total(3, i) = particle_set(i)%f(3) + f_nonbond(3, i)
END DO
IF (shell_present) THEN
ALLOCATE (fshell_total(3, nshell))
fshell_total = 0.0_dp
ALLOCATE (fcore_total(3, nshell))
fcore_total = 0.0_dp
DO i = 1, nshell
fcore_total(1, i) = core_particle_set(i)%f(1) + fcore_nonbond(1, i)
fcore_total(2, i) = core_particle_set(i)%f(2) + fcore_nonbond(2, i)
fcore_total(3, i) = core_particle_set(i)%f(3) + fcore_nonbond(3, i)
fshell_total(1, i) = shell_particle_set(i)%f(1) + fshell_nonbond(1, i)
fshell_total(2, i) = shell_particle_set(i)%f(2) + fshell_nonbond(2, i)
fshell_total(3, i) = shell_particle_set(i)%f(3) + fshell_nonbond(3, i)
END DO
END IF
IF (iw > 0) THEN
WRITE (iw, '(A)') " FIST::(1)INTERNAL + ELECTROSTATIC BONDED + NONBONDED"
WRITE (iw, '(3f15.9)') ((f_total(j, i), j=1, 3), i=1, natoms)
IF (shell_present .AND. shell_model_ad) THEN
WRITE (iw, '(A)') " FIST::(1)SHELL INTERNAL + ELECTROSTATIC BONDED + NONBONDED"
WRITE (iw, '(3f15.9)') ((fshell_total(j, i), j=1, 3), i=1, nshell)
WRITE (iw, '(A)') " FIST::(1)CORE INTERNAL + ELECTROSTATIC BONDED + NONBONDED"
WRITE (iw, '(3f15.9)') ((fcore_total(j, i), j=1, 3), i=1, nshell)
END IF
END IF
! Adding in the reciprocal forces: EWALD is a special case because of distrubted data
IF (ewald_type == do_ewald_ewald) THEN
node = 0
DO iparticle_kind = 1, nparticle_kind
nparticle_local = local_particles%n_el(iparticle_kind)
DO iparticle_local = 1, nparticle_local
ii = local_particles%list(iparticle_kind)%array(iparticle_local)
node = node + 1
f_total(1, ii) = f_total(1, ii) + fg_coulomb(1, node)
f_total(2, ii) = f_total(2, ii) + fg_coulomb(2, node)
f_total(3, ii) = f_total(3, ii) + fg_coulomb(3, node)
IF (PRESENT(debug)) THEN
debug%f_g(1, ii) = debug%f_g(1, ii) + fg_coulomb(1, node)
debug%f_g(2, ii) = debug%f_g(2, ii) + fg_coulomb(2, node)
debug%f_g(3, ii) = debug%f_g(3, ii) + fg_coulomb(3, node)
END IF
IF (atprop_env%energy) THEN
atprop_env%atener(ii) = atprop_env%atener(ii) + e_coulomb(node)
END IF
END DO
END DO
IF (atprop_env%energy) THEN
DEALLOCATE (e_coulomb)
END IF
END IF
IF (iw > 0) THEN
WRITE (iw, '(A)') " FIST::(2)TOTAL FORCES(1)+ ELECTROSTATIC FORCES"
WRITE (iw, '(3f15.9)') ((f_total(j, i), j=1, 3), i=1, natoms)
IF (shell_present .AND. shell_model_ad) THEN
WRITE (iw, '(A)') " FIST::(2)SHELL TOTAL FORCES(1)+ ELECTROSTATIC FORCES "
WRITE (iw, '(3f15.9)') ((fshell_total(j, i), j=1, 3), i=1, nshell)
WRITE (iw, '(A)') " FIST::(2)CORE TOTAL FORCES(1)+ ELECTROSTATIC FORCES"
WRITE (iw, '(3f15.9)') ((fcore_total(j, i), j=1, 3), i=1, nshell)
END IF
END IF
IF (use_virial) THEN
! Add up all the pressure tensors
IF (ewald_type == do_ewald_none) THEN
virial%pv_virial = pv_nonbond + pv_bond + pv_bend + pv_torsion + pv_imptors + pv_urey_bradley
CALL para_env%sum(virial%pv_virial)
ELSE
ident = 0.0_dp
DO i = 1, 3
ident(i, i) = 1.0_dp
END DO
virial%pv_virial = pv_nonbond + pv_bond + pv_bend + pv_torsion + pv_imptors + pv_urey_bradley + pv_bc
CALL para_env%sum(virial%pv_virial)
virial%pv_virial = virial%pv_virial + ident*thermo%e_neut
virial%pv_virial = virial%pv_virial + pv_g
END IF
END IF
! Sum total forces
CALL para_env%sum(f_total)
IF (shell_present .AND. shell_model_ad) THEN
CALL para_env%sum(fcore_total)
CALL para_env%sum(fshell_total)
END IF
! contributions from fields (currently all quantities are fully replicated!)
IF (efield%apply_field) THEN
thermo%pot = thermo%pot + ef_ener
f_total(1:3, 1:natoms) = f_total(1:3, 1:natoms) + ef_f(1:3, 1:natoms)
IF (use_virial) THEN
virial%pv_virial(1:3, 1:3) = virial%pv_virial(1:3, 1:3) + ef_pv(1:3, 1:3)
END IF
DEALLOCATE (ef_f)
END IF
! Assign to particle_set
SELECT CASE (ewald_type)
CASE (do_ewald_spme, do_ewald_pme)
IF (shell_present .AND. shell_model_ad) THEN
DO i = 1, natoms
shell_index = particle_set(i)%shell_index
IF (shell_index == 0) THEN
particle_set(i)%f(1:3) = f_total(1:3, i) + fg_coulomb(1:3, i)
ELSE
atomic_kind => particle_set(i)%atomic_kind
CALL get_atomic_kind(atomic_kind=atomic_kind, shell=shell, mass=mass)
fc = shell%mass_core/mass
fcore_total(1:3, shell_index) = fcore_total(1:3, shell_index) + particle_set(i)%f(1:3)*fc
fs = shell%mass_shell/mass
fshell_total(1:3, shell_index) = fshell_total(1:3, shell_index) + particle_set(i)%f(1:3)*fs
END IF
END DO
DO i = 1, nshell
core_particle_set(i)%f(1) = fcore_total(1, i) + fgcore_coulomb(1, i)
core_particle_set(i)%f(2) = fcore_total(2, i) + fgcore_coulomb(2, i)
core_particle_set(i)%f(3) = fcore_total(3, i) + fgcore_coulomb(3, i)
shell_particle_set(i)%f(1) = fshell_total(1, i) + fgshell_coulomb(1, i)
shell_particle_set(i)%f(2) = fshell_total(2, i) + fgshell_coulomb(2, i)
shell_particle_set(i)%f(3) = fshell_total(3, i) + fgshell_coulomb(3, i)
END DO
ELSEIF (shell_present .AND. .NOT. shell_model_ad) THEN
CPABORT("Non adiabatic shell-model not implemented.")
ELSE
DO i = 1, natoms
particle_set(i)%f(1) = f_total(1, i) + fg_coulomb(1, i)
particle_set(i)%f(2) = f_total(2, i) + fg_coulomb(2, i)
particle_set(i)%f(3) = f_total(3, i) + fg_coulomb(3, i)
END DO
END IF
CASE (do_ewald_ewald, do_ewald_none)
IF (shell_present .AND. shell_model_ad) THEN
DO i = 1, natoms
shell_index = particle_set(i)%shell_index
IF (shell_index == 0) THEN
particle_set(i)%f(1:3) = f_total(1:3, i)
ELSE
atomic_kind => particle_set(i)%atomic_kind
CALL get_atomic_kind(atomic_kind=atomic_kind, shell=shell, mass=mass)
fc = shell%mass_core/mass
fcore_total(1:3, shell_index) = fcore_total(1:3, shell_index) + particle_set(i)%f(1:3)*fc
fs = shell%mass_shell/mass
fshell_total(1:3, shell_index) = fshell_total(1:3, shell_index) + particle_set(i)%f(1:3)*fs
END IF
END DO
DO i = 1, nshell
core_particle_set(i)%f(1) = fcore_total(1, i)
core_particle_set(i)%f(2) = fcore_total(2, i)
core_particle_set(i)%f(3) = fcore_total(3, i)
shell_particle_set(i)%f(1) = fshell_total(1, i)
shell_particle_set(i)%f(2) = fshell_total(2, i)
shell_particle_set(i)%f(3) = fshell_total(3, i)
END DO
ELSEIF (shell_present .AND. .NOT. shell_model_ad) THEN
CPABORT("Non adiabatic shell-model not implemented.")
ELSE
DO i = 1, natoms
particle_set(i)%f(1) = f_total(1, i)
particle_set(i)%f(2) = f_total(2, i)
particle_set(i)%f(3) = f_total(3, i)
END DO
END IF
END SELECT
IF (iw > 0) THEN
WRITE (iw, '(A)') " FIST::(3)TOTAL FORCES - THE END..."
WRITE (iw, '(3f15.9)') ((particle_set(i)%f(j), j=1, 3), i=1, natoms)
IF (shell_present .AND. shell_model_ad) THEN
WRITE (iw, '(A)') " FIST::(3)SHELL TOTAL FORCES - THE END..."
WRITE (iw, '(3f15.9)') ((shell_particle_set(i)%f(j), j=1, 3), i=1, nshell)
WRITE (iw, '(A)') " FIST::(3)CORE TOTAL FORCES - THE END..."
WRITE (iw, '(3f15.9)') ((core_particle_set(i)%f(j), j=1, 3), i=1, nshell)
END IF
WRITE (iw, '(A,f15.9)') "Energy after FIST calculation.. exiting now ::", thermo%pot
END IF
!
! If we are doing debugging, check if variables are present and assign
!
IF (PRESENT(debug)) THEN
CALL para_env%sum(pot_manybody)
debug%pot_manybody = pot_manybody
CALL para_env%sum(pot_nonbond)
debug%pot_nonbond = pot_nonbond
CALL para_env%sum(pot_bond)
debug%pot_bond = pot_bond
CALL para_env%sum(pot_bend)
debug%pot_bend = pot_bend
CALL para_env%sum(pot_torsion)
debug%pot_torsion = pot_torsion
CALL para_env%sum(pot_imptors)
debug%pot_imptors = pot_imptors
CALL para_env%sum(pot_opbend)
debug%pot_opbend = pot_opbend
CALL para_env%sum(pot_urey_bradley)
debug%pot_urey_bradley = pot_urey_bradley
CALL para_env%sum(f_nonbond)
debug%f_nonbond = f_nonbond
IF (use_virial) THEN
CALL para_env%sum(pv_nonbond)
debug%pv_nonbond = pv_nonbond
CALL para_env%sum(pv_bond)
debug%pv_bond = pv_bond
CALL para_env%sum(pv_bend)
debug%pv_bend = pv_bend
CALL para_env%sum(pv_torsion)
debug%pv_torsion = pv_torsion
CALL para_env%sum(pv_imptors)
debug%pv_imptors = pv_imptors
CALL para_env%sum(pv_urey_bradley)
debug%pv_ub = pv_urey_bradley
END IF
SELECT CASE (ewald_type)
CASE (do_ewald_spme, do_ewald_pme)
debug%pot_g = vg_coulomb
debug%pv_g = pv_g
debug%f_g = fg_coulomb
CASE (do_ewald_ewald)
debug%pot_g = vg_coulomb
IF (use_virial) debug%pv_g = pv_g
CASE default
debug%pot_g = 0.0_dp
debug%f_g = 0.0_dp
IF (use_virial) debug%pv_g = 0.0_dp
END SELECT
END IF
! print properties if requested
print_section => section_vals_get_subs_vals(mm_section, "PRINT")
CALL print_fist(fist_env, print_section, atomic_kind_set, particle_set, cell)
! deallocating all local variables
IF (ALLOCATED(fg_coulomb)) THEN
DEALLOCATE (fg_coulomb)
END IF
IF (ALLOCATED(f_total)) THEN
DEALLOCATE (f_total)
END IF
DEALLOCATE (f_nonbond)
IF (shell_present) THEN
DEALLOCATE (fshell_total)
IF (ewald_type /= do_ewald_none) THEN
DEALLOCATE (fgshell_coulomb)
END IF
DEALLOCATE (fshell_nonbond)
END IF
CALL cp_print_key_finished_output(iw, logger, mm_section, "PRINT%DERIVATIVES")
CALL cp_print_key_finished_output(iw2, logger, mm_section, "PRINT%EWALD_INFO")
CALL timestop(handle)
END SUBROUTINE fist_calc_energy_force
! **************************************************************************************************
!> \brief Print properties number according the requests in input file
!> \param fist_env ...
!> \param print_section ...
!> \param atomic_kind_set ...
!> \param particle_set ...
!> \param cell ...
!> \par History
!> [01.2006] created
!> \author Teodoro Laino
! **************************************************************************************************
SUBROUTINE print_fist(fist_env, print_section, atomic_kind_set, particle_set, cell)
TYPE(fist_environment_type), POINTER :: fist_env
TYPE(section_vals_type), POINTER :: print_section
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(cell_type), POINTER :: cell
INTEGER :: unit_nr
TYPE(cp_logger_type), POINTER :: logger
TYPE(fist_nonbond_env_type), POINTER :: fist_nonbond_env
TYPE(section_vals_type), POINTER :: print_key
NULLIFY (logger, print_key, fist_nonbond_env)
logger => cp_get_default_logger()
print_key => section_vals_get_subs_vals(print_section, "dipole")
CALL fist_env_get(fist_env, fist_nonbond_env=fist_nonbond_env)
IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), &
cp_p_file)) THEN
unit_nr = cp_print_key_unit_nr(logger, print_section, "dipole", &
extension=".data", middle_name="MM_DIPOLE", log_filename=.FALSE.)
CALL fist_dipole(fist_env, print_section, atomic_kind_set, particle_set, &
cell, unit_nr, fist_nonbond_env%charges)
CALL cp_print_key_finished_output(unit_nr, logger, print_key)
END IF
END SUBROUTINE print_fist
END MODULE fist_force