LCOV - code coverage report
Current view: top level - src - manybody_gal21.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 80.1 % 332 266
Test Date: 2025-07-25 12:55:17 Functions: 90.0 % 10 9

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : ! **************************************************************************************************
       9              : !> \brief Implementation of the GAL21 potential
      10              : !>
      11              : !> \author Clabaut Paul
      12              : ! **************************************************************************************************
      13              : MODULE manybody_gal21
      14              : 
      15              :    USE atomic_kind_types,               ONLY: get_atomic_kind
      16              :    USE cell_types,                      ONLY: cell_type,&
      17              :                                               pbc
      18              :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      19              :                                               cp_logger_type,&
      20              :                                               cp_to_string
      21              :    USE cp_output_handling,              ONLY: cp_print_key_finished_output,&
      22              :                                               cp_print_key_unit_nr
      23              :    USE fist_neighbor_list_types,        ONLY: fist_neighbor_type,&
      24              :                                               neighbor_kind_pairs_type
      25              :    USE fist_nonbond_env_types,          ONLY: pos_type
      26              :    USE input_section_types,             ONLY: section_vals_type
      27              :    USE kinds,                           ONLY: dp
      28              :    USE message_passing,                 ONLY: mp_para_env_type
      29              :    USE pair_potential_types,            ONLY: gal21_pot_type,&
      30              :                                               gal21_type,&
      31              :                                               pair_potential_pp_type,&
      32              :                                               pair_potential_single_type
      33              :    USE particle_types,                  ONLY: particle_type
      34              :    USE util,                            ONLY: sort
      35              : #include "./base/base_uses.f90"
      36              : 
      37              :    IMPLICIT NONE
      38              : 
      39              :    PRIVATE
      40              :    PUBLIC :: setup_gal21_arrays, destroy_gal21_arrays, &
      41              :              gal21_energy, gal21_forces, &
      42              :              print_nr_ions_gal21
      43              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'manybody_gal21'
      44              : 
      45              : CONTAINS
      46              : 
      47              : ! **************************************************************************************************
      48              : !> \brief  Main part of the energy evaluation of GAL2119
      49              : !> \param pot_loc value of total potential energy
      50              : !> \param gal21 all parameters of GAL2119
      51              : !> \param r_last_update_pbc position of every atoms on previous frame
      52              : !> \param iparticle first index of the atom of the evaluated pair
      53              : !> \param jparticle second index of the atom of the evaluated pair
      54              : !> \param cell dimension of the pbc cell
      55              : !> \param particle_set full list of atoms of the system
      56              : !> \param mm_section ...
      57              : !> \author Clabaut Paul - 2019 - ENS de Lyon
      58              : ! **************************************************************************************************
      59         5732 :    SUBROUTINE gal21_energy(pot_loc, gal21, r_last_update_pbc, iparticle, jparticle, &
      60              :                            cell, particle_set, mm_section)
      61              : 
      62              :       REAL(KIND=dp), INTENT(OUT)                         :: pot_loc
      63              :       TYPE(gal21_pot_type), POINTER                      :: gal21
      64              :       TYPE(pos_type), DIMENSION(:), POINTER              :: r_last_update_pbc
      65              :       INTEGER, INTENT(IN)                                :: iparticle, jparticle
      66              :       TYPE(cell_type), POINTER                           :: cell
      67              :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
      68              :       TYPE(section_vals_type), POINTER                   :: mm_section
      69              : 
      70              :       CHARACTER(LEN=2)                                   :: element_symbol
      71              :       INTEGER                                            :: index_outfile
      72              :       REAL(KIND=dp)                                      :: anglepart, AO, BO, bxy, bz, cosalpha, &
      73              :                                                             drji2, eps, nvec(3), rji(3), sinalpha, &
      74              :                                                             sum_weight, Vang, Vgaussian, VH, VTT, &
      75              :                                                             weight
      76              :       TYPE(cp_logger_type), POINTER                      :: logger
      77              : 
      78         5732 :       pot_loc = 0.0_dp
      79              :       CALL get_atomic_kind(atomic_kind=particle_set(iparticle)%atomic_kind, &
      80         5732 :                            element_symbol=element_symbol) !Read the atom type of i
      81              : 
      82         5732 :       IF (element_symbol == "O") THEN !To avoid counting two times each pair
      83              : 
      84         2866 :          rji(:) = pbc(r_last_update_pbc(jparticle)%r(:), r_last_update_pbc(iparticle)%r(:), cell) !Vector in pbc from j to i
      85              : 
      86         2866 :          IF (.NOT. ALLOCATED(gal21%n_vectors)) THEN !First calling of the forcefield only
      87            3 :             ALLOCATE (gal21%n_vectors(3, SIZE(particle_set)))
      88         3481 :             gal21%n_vectors(:, :) = 0.0_dp
      89              :          END IF
      90              : 
      91         2866 :          IF (gal21%express) THEN
      92            0 :             logger => cp_get_default_logger()
      93              :             index_outfile = cp_print_key_unit_nr(logger, mm_section, &
      94            0 :                                                  "PRINT%PROGRAM_RUN_INFO", extension=".mmLog")
      95            0 :             IF (index_outfile > 0) WRITE (index_outfile, *) "GCN", gal21%gcn(jparticle)
      96              :             CALL cp_print_key_finished_output(index_outfile, logger, mm_section, &
      97            0 :                                               "PRINT%PROGRAM_RUN_INFO")
      98              :          END IF
      99              : 
     100              :          !Build epsilon attraction and the parameters of the gaussian attraction as a function of gcn
     101         2866 :          eps = gal21%epsilon1*gal21%gcn(jparticle)*gal21%gcn(jparticle) + gal21%epsilon2*gal21%gcn(jparticle) + gal21%epsilon3
     102         2866 :          bxy = gal21%bxy1 + gal21%bxy2*gal21%gcn(jparticle)
     103         2866 :          bz = gal21%bz1 + gal21%bz2*gal21%gcn(jparticle)
     104              : 
     105              :          !Angular dependance %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     106         2866 :          Vang = 0.0_dp
     107              : 
     108              :          !Calculation of the normal vector centered on the Me atom of the pair, only the first time that an interaction with the metal atom of the pair is evaluated
     109              :          IF (gal21%n_vectors(1, jparticle) == 0.0_dp .AND. &
     110         2866 :              gal21%n_vectors(2, jparticle) == 0.0_dp .AND. &
     111              :              gal21%n_vectors(3, jparticle) == 0.0_dp) THEN
     112              :             gal21%n_vectors(:, jparticle) = normale(gal21, r_last_update_pbc, jparticle, &
     113          596 :                                                     particle_set, cell)
     114              :          END IF
     115              : 
     116        11464 :          nvec(:) = gal21%n_vectors(:, jparticle) !Else, retrive it, should not have moved sinc metal is supposed to be frozen
     117              : 
     118              :          !Calculation of the sum of the expontial weights of each Me surrounding the principal one
     119         2866 :          sum_weight = somme(gal21, r_last_update_pbc, iparticle, particle_set, cell)
     120              : 
     121              :          !Exponential damping weight for angular dependance
     122        11464 :          weight = EXP(-SQRT(DOT_PRODUCT(rji, rji))/gal21%r1)
     123              : 
     124              :          !Calculation of the truncated fourier series of the water-dipole/surface-normal angle
     125              :          anglepart = 0.0_dp
     126              :          VH = 0.0_dp
     127              :          CALL angular(anglepart, VH, gal21, r_last_update_pbc, iparticle, jparticle, cell, particle_set, nvec, &
     128         2866 :                       .TRUE., mm_section)
     129              : 
     130              :          !Build the complete angular potential while avoiding division by 0
     131         2866 :          IF (weight /= 0) THEN
     132         2866 :             Vang = weight*weight*anglepart/sum_weight
     133         2866 :             IF (gal21%express) THEN
     134            0 :                logger => cp_get_default_logger()
     135              :                index_outfile = cp_print_key_unit_nr(logger, mm_section, &
     136            0 :                                                     "PRINT%PROGRAM_RUN_INFO", extension=".mmLog")
     137            0 :                IF (index_outfile > 0) WRITE (index_outfile, *) "Fermi", weight*weight/sum_weight
     138              :                CALL cp_print_key_finished_output(index_outfile, logger, mm_section, &
     139            0 :                                                  "PRINT%PROGRAM_RUN_INFO")
     140              :             END IF
     141              :          END IF
     142              :          !END Angular%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     143              : 
     144              :          !Attractive Gaussian %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     145         2866 :          Vgaussian = 0.0_dp
     146        11464 :          drji2 = DOT_PRODUCT(rji, rji)
     147              :          !Alpha is the angle of the Me-O vector with the normale vector. Used for gaussian attaction
     148              : 
     149        11464 :          cosalpha = DOT_PRODUCT(rji, nvec)/SQRT(drji2)
     150         2866 :          IF (cosalpha < -1.0_dp) cosalpha = -1.0_dp
     151         2866 :          IF (cosalpha > +1.0_dp) cosalpha = +1.0_dp
     152         2866 :          sinalpha = SIN(ACOS(cosalpha))
     153              : 
     154              :          !Gaussian component of the energy
     155              :          Vgaussian = -1.0_dp*eps*EXP(-bz*drji2*cosalpha*cosalpha &
     156         2866 :                                      - bxy*drji2*sinalpha*sinalpha)
     157              :          !END Gaussian%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     158              : 
     159         2866 :          AO = gal21%AO1 + gal21%AO2*gal21%gcn(jparticle)
     160         2866 :          BO = gal21%BO1 + gal21%BO2*gal21%gcn(jparticle)
     161              : 
     162              :          !Tang and toennies potential for physisorption
     163              :          VTT = AO*EXP(-BO*SQRT(drji2)) - (1.0 - EXP(-BO*SQRT(drji2)) &
     164              :                                           - BO*SQRT(drji2)*EXP(-BO*SQRT(drji2)) &
     165              :                                           - (((BO*SQRT(drji2))**2)/2)*EXP(-BO*SQRT(drji2)) &
     166              :                                           - (((BO*SQRT(drji2))**3)/6)*EXP(-BO*SQRT(drji2)) &
     167              :                                           - (((BO*SQRT(drji2))**4)/24)*EXP(-BO*SQRT(drji2)) &
     168              :                                           - (((BO*SQRT(drji2))**5)/120)*EXP(-BO*SQRT(drji2)) &
     169              :                                           - (((BO*SQRT(drji2))**6)/720)*EXP(-BO*SQRT(drji2))) &
     170         2866 :                *gal21%c/(SQRT(drji2)**6)
     171              : 
     172              :          !For fit purpose only
     173         2866 :          IF (gal21%express) THEN
     174            0 :             logger => cp_get_default_logger()
     175              :             index_outfile = cp_print_key_unit_nr(logger, mm_section, &
     176            0 :                                                  "PRINT%PROGRAM_RUN_INFO", extension=".mmLog")
     177            0 :             IF (index_outfile > 0) WRITE (index_outfile, *) "Gau", -1.0_dp*EXP(-bz*drji2*cosalpha*cosalpha &
     178            0 :                                                                                - bxy*drji2*sinalpha*sinalpha)
     179            0 :             IF (weight == 0 .AND. index_outfile > 0) WRITE (index_outfile, *) "Fermi  0"
     180            0 :             IF (index_outfile > 0) WRITE (index_outfile, *) "expO", EXP(-BO*SQRT(drji2))
     181            0 :             IF (index_outfile > 0) WRITE (index_outfile, *) "cstpart", -(1.0 - EXP(-BO*SQRT(drji2)) &
     182              :                                                                          - BO*SQRT(drji2)*EXP(-BO*SQRT(drji2)) &
     183              :                                                                          - (((BO*SQRT(drji2))**2)/2)*EXP(-BO*SQRT(drji2)) &
     184              :                                                                          - (((BO*SQRT(drji2))**3)/6)*EXP(-BO*SQRT(drji2)) &
     185              :                                                                          - (((BO*SQRT(drji2))**4)/24)*EXP(-BO*SQRT(drji2)) &
     186              :                                                                          - (((BO*SQRT(drji2))**5)/120)*EXP(-BO*SQRT(drji2)) &
     187              :                                                                          - (((BO*SQRT(drji2))**6)/720)*EXP(-BO*SQRT(drji2))) &
     188            0 :                *gal21%c/(SQRT(drji2)**6)
     189            0 :             IF (index_outfile > 0) WRITE (index_outfile, *) "params_lin_eps", gal21%epsilon1, gal21%epsilon2, gal21%epsilon3
     190            0 :             IF (index_outfile > 0) WRITE (index_outfile, *) "params_lin_A0", AO
     191              :             CALL cp_print_key_finished_output(index_outfile, logger, mm_section, &
     192            0 :                                               "PRINT%PROGRAM_RUN_INFO")
     193              :          END IF
     194              :          !Compute the total energy
     195         2866 :          pot_loc = Vgaussian + Vang + VTT + VH
     196              : 
     197              :       END IF
     198              : 
     199         5732 :    END SUBROUTINE gal21_energy
     200              : 
     201              : ! **************************************************************************************************
     202              : !> \brief The idea is to build a vector normal to the local surface by using the symetry of the
     203              : !>        surface that make the opposite vectors compensate themself. The vector is therefore in the
     204              : !>.       direction of the missing atoms of a large coordination sphere
     205              : !> \param gal21 ...
     206              : !> \param r_last_update_pbc ...
     207              : !> \param jparticle ...
     208              : !> \param particle_set ...
     209              : !> \param cell ...
     210              : !> \return ...
     211              : !> \retval normale ...
     212              : ! **************************************************************************************************
     213          149 :    FUNCTION normale(gal21, r_last_update_pbc, jparticle, particle_set, cell)
     214              :       TYPE(gal21_pot_type), POINTER                      :: gal21
     215              :       TYPE(pos_type), DIMENSION(:), POINTER              :: r_last_update_pbc
     216              :       INTEGER, INTENT(IN)                                :: jparticle
     217              :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     218              :       TYPE(cell_type), POINTER                           :: cell
     219              :       REAL(KIND=dp)                                      :: normale(3)
     220              : 
     221              :       CHARACTER(LEN=2)                                   :: element_symbol_k
     222              :       INTEGER                                            :: kparticle, natom
     223              :       REAL(KIND=dp)                                      :: drjk, rjk(3)
     224              : 
     225          149 :       natom = SIZE(particle_set)
     226          596 :       normale(:) = 0.0_dp
     227              : 
     228       129779 :       DO kparticle = 1, natom !Loop on every atom of the system
     229       129630 :          IF (kparticle == jparticle) CYCLE !Avoid the principal Me atom (j) in the counting
     230              :          CALL get_atomic_kind(atomic_kind=particle_set(kparticle)%atomic_kind, &
     231       129481 :                               element_symbol=element_symbol_k)
     232       129481 :          IF (element_symbol_k /= gal21%met1 .AND. element_symbol_k /= gal21%met2) CYCLE !Keep only metals
     233        28459 :          rjk(:) = pbc(r_last_update_pbc(jparticle)%r(:), r_last_update_pbc(kparticle)%r(:), cell)
     234       113836 :          drjk = SQRT(DOT_PRODUCT(rjk, rjk))
     235        28459 :          IF (drjk > gal21%rcutsq) CYCLE !Keep only those within square root of the force-field cutoff distance of the metallic atom of the evaluated pair
     236       215156 :          normale(:) = normale(:) - rjk(:)/(drjk*drjk*drjk*drjk*drjk) !Build the normal, vector by vector
     237              :       END DO
     238              : 
     239              :       ! Normalisation of the vector
     240         1043 :       normale(:) = normale(:)/SQRT(DOT_PRODUCT(normale, normale))
     241              : 
     242              :    END FUNCTION normale
     243              : 
     244              : ! **************************************************************************************************
     245              : !> \brief Scan all the Me atoms that have been counted in the O-Me paires and sum their exp. weights
     246              : !> \param gal21 ...
     247              : !> \param r_last_update_pbc ...
     248              : !> \param iparticle ...
     249              : !> \param particle_set ...
     250              : !> \param cell ...
     251              : !> \return ...
     252              : !> \retval somme ...
     253              : ! **************************************************************************************************
     254         5732 :    FUNCTION somme(gal21, r_last_update_pbc, iparticle, particle_set, cell)
     255              :       TYPE(gal21_pot_type), POINTER                      :: gal21
     256              :       TYPE(pos_type), DIMENSION(:), POINTER              :: r_last_update_pbc
     257              :       INTEGER, INTENT(IN)                                :: iparticle
     258              :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     259              :       TYPE(cell_type), POINTER                           :: cell
     260              :       REAL(KIND=dp)                                      :: somme
     261              : 
     262              :       CHARACTER(LEN=2)                                   :: element_symbol_k
     263              :       INTEGER                                            :: kparticle, natom
     264              :       REAL(KIND=dp)                                      :: rki(3)
     265              : 
     266         5732 :       natom = SIZE(particle_set)
     267         5732 :       somme = 0.0_dp
     268              : 
     269      4992572 :       DO kparticle = 1, natom !Loop on every atom of the system
     270              :          CALL get_atomic_kind(atomic_kind=particle_set(kparticle)%atomic_kind, &
     271      4986840 :                               element_symbol=element_symbol_k)
     272      4986840 :          IF (element_symbol_k /= gal21%met1 .AND. element_symbol_k /= gal21%met2) CYCLE !Keep only metals
     273      1100544 :          rki(:) = pbc(r_last_update_pbc(kparticle)%r(:), r_last_update_pbc(iparticle)%r(:), cell)
     274      4402176 :          IF (SQRT(DOT_PRODUCT(rki, rki)) > gal21%rcutsq) CYCLE !Keep only those within cutoff distance of the oxygen atom of the evaluated pair (the omega ensemble)
     275      4402176 :          IF (element_symbol_k == gal21%met1) somme = somme + EXP(-SQRT(DOT_PRODUCT(rki, rki))/gal21%r1) !Build the sum of the exponential weights
     276      1106276 :          IF (element_symbol_k == gal21%met2) somme = somme + EXP(-SQRT(DOT_PRODUCT(rki, rki))/gal21%r2) !Build the sum of the exponential weights
     277              :       END DO
     278              : 
     279         5732 :    END FUNCTION somme
     280              : 
     281              : ! **************************************************************************************************
     282              : !> \brief Compute the angular dependance (on theta) of the forcefield
     283              : !> \param anglepart ...
     284              : !> \param VH ...
     285              : !> \param gal21 ...
     286              : !> \param r_last_update_pbc ...
     287              : !> \param iparticle ...
     288              : !> \param jparticle ...
     289              : !> \param cell ...
     290              : !> \param particle_set ...
     291              : !> \param nvec ...
     292              : !> \param energy ...
     293              : !> \param mm_section ...
     294              : !> \return ...
     295              : !> \retval angular ...
     296              : ! **************************************************************************************************
     297         5732 :    SUBROUTINE angular(anglepart, VH, gal21, r_last_update_pbc, iparticle, jparticle, cell, &
     298              :                       particle_set, nvec, energy, mm_section)
     299              :       REAL(KIND=dp)                                      :: anglepart, VH
     300              :       TYPE(gal21_pot_type), POINTER                      :: gal21
     301              :       TYPE(pos_type), DIMENSION(:), POINTER              :: r_last_update_pbc
     302              :       INTEGER, INTENT(IN)                                :: iparticle, jparticle
     303              :       TYPE(cell_type), POINTER                           :: cell
     304              :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     305              :       REAL(KIND=dp), DIMENSION(3)                        :: nvec
     306              :       LOGICAL                                            :: energy
     307              :       TYPE(section_vals_type), POINTER                   :: mm_section
     308              : 
     309              :       CHARACTER(LEN=2)                                   :: element_symbol
     310              :       INTEGER                                            :: count_h, iatom, index_h1, index_h2, &
     311              :                                                             index_outfile, natom
     312              :       REAL(KIND=dp)                                      :: a1, a2, a3, a4, BH, costheta, &
     313              :                                                             h_max_dist, rih(3), rih1(3), rih2(3), &
     314              :                                                             rix(3), rjh1(3), rjh2(3), theta
     315              :       TYPE(cp_logger_type), POINTER                      :: logger
     316              : 
     317         5732 :       count_h = 0
     318         5732 :       index_h1 = 0
     319         5732 :       index_h2 = 0
     320         5732 :       h_max_dist = 2.1_dp ! 1.1 angstrom
     321         5732 :       natom = SIZE(particle_set)
     322              : 
     323      4992572 :       DO iatom = 1, natom !Loop on every atom of the system
     324              :          CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
     325      4986840 :                               element_symbol=element_symbol)
     326      4986840 :          IF (element_symbol /= "H") CYCLE !Kepp only hydrogen
     327      2590864 :          rih(:) = pbc(r_last_update_pbc(iparticle)%r(:), r_last_update_pbc(iatom)%r(:), cell)
     328     10363456 :          IF (SQRT(DOT_PRODUCT(rih, rih)) >= h_max_dist) CYCLE !Keep only hydrogen that are bounded to the considered O
     329        11464 :          count_h = count_h + 1
     330        17196 :          IF (count_h == 1) THEN
     331              :             index_h1 = iatom
     332         5732 :          ELSEIF (count_h == 2) THEN
     333         5732 :             index_h2 = iatom
     334              :          END IF
     335              :       END DO
     336              : 
     337              :       ! Abort if the oxygen is not part of a water molecule (2 H)
     338         5732 :       IF (count_h /= 2) THEN
     339              :          CALL cp_abort(__LOCATION__, &
     340            0 :                        " Error: Found "//cp_to_string(count_h)//" H atoms for O atom "//cp_to_string(iparticle))
     341              :       END IF
     342              : 
     343         5732 :       a1 = gal21%a11 + gal21%a12*gal21%gcn(jparticle) + gal21%a13*gal21%gcn(jparticle)*gal21%gcn(jparticle)
     344         5732 :       a2 = gal21%a21 + gal21%a22*gal21%gcn(jparticle) + gal21%a23*gal21%gcn(jparticle)*gal21%gcn(jparticle)
     345         5732 :       a3 = gal21%a31 + gal21%a32*gal21%gcn(jparticle) + gal21%a33*gal21%gcn(jparticle)*gal21%gcn(jparticle)
     346         5732 :       a4 = gal21%a41 + gal21%a42*gal21%gcn(jparticle) + gal21%a43*gal21%gcn(jparticle)*gal21%gcn(jparticle)
     347              : 
     348         5732 :       rih1(:) = pbc(r_last_update_pbc(iparticle)%r(:), r_last_update_pbc(index_h1)%r(:), cell)
     349         5732 :       rih2(:) = pbc(r_last_update_pbc(iparticle)%r(:), r_last_update_pbc(index_h2)%r(:), cell)
     350        22928 :       rix(:) = rih1(:) + rih2(:) ! build the dipole vector rix of the H2O molecule
     351        40124 :       costheta = DOT_PRODUCT(rix, nvec)/SQRT(DOT_PRODUCT(rix, rix))
     352         5732 :       IF (costheta < -1.0_dp) costheta = -1.0_dp
     353         5732 :       IF (costheta > +1.0_dp) costheta = +1.0_dp
     354         5732 :       theta = ACOS(costheta) ! Theta is the angle between the normal to the surface and the dipole
     355              :       anglepart = a1*costheta + a2*COS(2.0_dp*theta) + a3*COS(3.0_dp*theta) &
     356         5732 :                   + a4*COS(4.0_dp*theta) ! build the fourier series
     357              : 
     358         5732 :       BH = gal21%BH1 + gal21%gcn(jparticle)*gal21%BH2
     359              : 
     360         5732 :       rjh1(:) = pbc(r_last_update_pbc(jparticle)%r(:), r_last_update_pbc(index_h1)%r(:), cell)
     361         5732 :       rjh2(:) = pbc(r_last_update_pbc(jparticle)%r(:), r_last_update_pbc(index_h2)%r(:), cell)
     362              :       VH = (gal21%AH2*gal21%gcn(jparticle) + gal21%AH1)*(EXP(-BH*SQRT(DOT_PRODUCT(rjh1, rjh1))) + &
     363        40124 :                                                          EXP(-BH*SQRT(DOT_PRODUCT(rjh2, rjh2))))
     364              : 
     365              :       ! For fit purpose
     366         5732 :       IF (gal21%express .AND. energy) THEN
     367            0 :          logger => cp_get_default_logger()
     368              :          index_outfile = cp_print_key_unit_nr(logger, mm_section, &
     369            0 :                                               "PRINT%PROGRAM_RUN_INFO", extension=".mmLog")
     370              : 
     371            0 :          IF (index_outfile > 0) WRITE (index_outfile, *) "Fourier", costheta, COS(2.0_dp*theta), COS(3.0_dp*theta), &
     372            0 :             COS(4.0_dp*theta) !, theta
     373            0 :          IF (index_outfile > 0) WRITE (index_outfile, *) "H_rep", EXP(-BH*SQRT(DOT_PRODUCT(rjh1, rjh1))) + &
     374            0 :             EXP(-BH*SQRT(DOT_PRODUCT(rjh2, rjh2)))
     375              :          !IF (index_outfile > 0) WRITE (index_outfile, *) "H_r6", -1/DOT_PRODUCT(rjh1,rjh1)**3 -1/DOT_PRODUCT(rjh2,rjh2)**3
     376              : 
     377              :          CALL cp_print_key_finished_output(index_outfile, logger, mm_section, &
     378            0 :                                            "PRINT%PROGRAM_RUN_INFO")
     379              :       END IF
     380              : 
     381         5732 :    END SUBROUTINE
     382              : 
     383              : ! **************************************************************************************************
     384              : !> \brief forces generated by the GAL2119 potential
     385              : !> \param gal21 all parameters of GAL2119
     386              : !> \param r_last_update_pbc position of every atoms on previous frame
     387              : !> \param iparticle first index of the atom of the evaluated pair
     388              : !> \param jparticle second index of the atom of the evaluated pair
     389              : !> \param f_nonbond all the forces applying on the system
     390              : !> \param pv_nonbond ...
     391              : !> \param use_virial request of usage of virial (for barostat)
     392              : !> \param cell dimension of the pbc cell
     393              : !> \param particle_set full list of atoms of the system
     394              : !> \author Clabaut Paul - 2019 - ENS de Lyon
     395              : ! **************************************************************************************************
     396         5732 :    SUBROUTINE gal21_forces(gal21, r_last_update_pbc, iparticle, jparticle, f_nonbond, pv_nonbond, &
     397              :                            use_virial, cell, particle_set)
     398              :       TYPE(gal21_pot_type), POINTER                      :: gal21
     399              :       TYPE(pos_type), DIMENSION(:), POINTER              :: r_last_update_pbc
     400              :       INTEGER, INTENT(IN)                                :: iparticle, jparticle
     401              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT)      :: f_nonbond, pv_nonbond
     402              :       LOGICAL, INTENT(IN)                                :: use_virial
     403              :       TYPE(cell_type), POINTER                           :: cell
     404              :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     405              : 
     406              :       CHARACTER(LEN=2)                                   :: element_symbol
     407              :       REAL(KIND=dp) :: anglepart, AO, BO, bxy, bz, cosalpha, dGauss(3), drji, drjicosalpha(3), &
     408              :          drjisinalpha(3), dTT(3), dweight(3), eps, nvec(3), prefactor, rji(3), rji_hat(3), &
     409              :          sinalpha, sum_weight, Vgaussian, VH, weight
     410              :       TYPE(section_vals_type), POINTER                   :: mm_section
     411              : 
     412              :       CALL get_atomic_kind(atomic_kind=particle_set(iparticle)%atomic_kind, &
     413         5732 :                            element_symbol=element_symbol)
     414              : 
     415         5732 :       IF (element_symbol == "O") THEN !To avoid counting two times each pair
     416              : 
     417         2866 :          rji(:) = pbc(r_last_update_pbc(jparticle)%r(:), r_last_update_pbc(iparticle)%r(:), cell)
     418        11464 :          drji = SQRT(DOT_PRODUCT(rji, rji))
     419        11464 :          rji_hat(:) = rji(:)/drji ! hat = pure directional component of a given vector
     420              : 
     421         2866 :          IF (.NOT. ALLOCATED(gal21%n_vectors)) THEN !First calling of the forcefield only
     422            0 :             ALLOCATE (gal21%n_vectors(3, SIZE(particle_set)))
     423            0 :             gal21%n_vectors(:, :) = 0.0_dp
     424              :          END IF
     425              : 
     426              :          !Build epsilon attraction and the a parameters of the Fourier serie as quadratic fucntion of gcn
     427         2866 :          eps = gal21%epsilon1*gal21%gcn(jparticle)*gal21%gcn(jparticle) + gal21%epsilon2*gal21%gcn(jparticle) + gal21%epsilon3
     428         2866 :          bxy = gal21%bxy1 + gal21%bxy2*gal21%gcn(jparticle)
     429         2866 :          bz = gal21%bz1 + gal21%bz2*gal21%gcn(jparticle)
     430              : 
     431              :          !Angular dependance %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     432              : 
     433              :          !Calculation of the normal vector centered on the Me atom of the pair, only the first time that an interaction with the metal atom of the pair is evaluated
     434              :          IF (gal21%n_vectors(1, jparticle) == 0.0_dp .AND. &
     435         2866 :              gal21%n_vectors(2, jparticle) == 0.0_dp .AND. &
     436              :              gal21%n_vectors(3, jparticle) == 0.0_dp) THEN
     437              :             gal21%n_vectors(:, jparticle) = normale(gal21, r_last_update_pbc, jparticle, &
     438            0 :                                                     particle_set, cell)
     439              :          END IF
     440              : 
     441        11464 :          nvec(:) = gal21%n_vectors(:, jparticle) !Else, retrive it, should not have moved sinc metal is supposed to be frozen
     442              : 
     443              :          !Calculation of the sum of the expontial weights of each Me surrounding the principal one
     444         2866 :          sum_weight = somme(gal21, r_last_update_pbc, iparticle, particle_set, cell)
     445              : 
     446              :          !Exponential damping weight for angular dependance
     447         2866 :          weight = EXP(-drji/gal21%r1)
     448        11464 :          dweight(:) = 1.0_dp/gal21%r1*weight*rji_hat(:) !Derivativ of it
     449              : 
     450              :          !Calculation of the truncated fourier series of the water-dipole/surface-normal angle
     451         2866 :          NULLIFY (mm_section)
     452              :          anglepart = 0.0_dp
     453              :          VH = 0.0_dp
     454              :          CALL angular(anglepart, VH, gal21, r_last_update_pbc, iparticle, jparticle, cell, particle_set, nvec, &
     455         2866 :                       .FALSE., mm_section)
     456              : 
     457              :          !Build the average of the exponential weight while avoiding division by 0
     458         2866 :          IF (weight /= 0) THEN
     459              :             ! Calculate the first component of the derivativ of the angular term
     460              :             f_nonbond(1:3, iparticle) = f_nonbond(1:3, iparticle) + 2.0_dp*dweight(1:3)*weight* &
     461        11464 :                                         anglepart/sum_weight
     462              : 
     463         2866 :             IF (use_virial) THEN
     464              :                pv_nonbond(1, 1:3) = pv_nonbond(1, 1:3) + rji(1)*2.0_dp*dweight(1:3)*weight* &
     465            0 :                                     anglepart/sum_weight
     466              :                pv_nonbond(2, 1:3) = pv_nonbond(2, 1:3) + rji(2)*2.0_dp*dweight(1:3)*weight* &
     467            0 :                                     anglepart/sum_weight
     468              :                pv_nonbond(3, 1:3) = pv_nonbond(3, 1:3) + rji(3)*2.0_dp*dweight(1:3)*weight* &
     469            0 :                                     anglepart/sum_weight
     470              :             END IF
     471              : 
     472              :             ! Calculate the second component of the derivativ of the angular term
     473              :             CALL somme_d(gal21, r_last_update_pbc, iparticle, jparticle, &
     474         2866 :                          f_nonbond, pv_nonbond, use_virial, particle_set, cell, anglepart, sum_weight)
     475              : 
     476         2866 :             prefactor = (-1.0_dp)*weight*weight/sum_weight ! Avoiding division by 0
     477              : 
     478              :             ! Calculate the third component of the derivativ of the angular term
     479              :             CALL angular_d(gal21, r_last_update_pbc, iparticle, jparticle, &
     480         2866 :                            f_nonbond, pv_nonbond, use_virial, prefactor, cell, particle_set, nvec)
     481              :          END IF
     482              :          !END Angular%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     483              : 
     484              :          !Attractive Gaussian %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     485              :          !Alpha is the angle of the Me-O vector with the normale vector. Used for gaussian attaction
     486        11464 :          cosalpha = DOT_PRODUCT(rji, nvec)/drji
     487         2866 :          IF (cosalpha < -1.0_dp) cosalpha = -1.0_dp
     488         2866 :          IF (cosalpha > +1.0_dp) cosalpha = +1.0_dp
     489         2866 :          sinalpha = SIN(ACOS(cosalpha))
     490              : 
     491              :          !Gaussian component of the energy
     492              :          Vgaussian = -1.0_dp*eps*EXP(-bz*DOT_PRODUCT(rji, rji)*cosalpha*cosalpha &
     493        11464 :                                      - bxy*DOT_PRODUCT(rji, rji)*sinalpha*sinalpha)
     494              : 
     495              :          ! Calculation of partial derivativ of the gaussian components
     496        11464 :          drjicosalpha(:) = rji_hat(:)*cosalpha + nvec(:) - cosalpha*rji_hat(:)
     497        11464 :          drjisinalpha(:) = rji_hat(:)*sinalpha - (cosalpha/sinalpha)*(nvec(:) - cosalpha*rji_hat(:))
     498              :          dGauss(:) = (-1.0_dp*bz*2*drji*cosalpha*drjicosalpha - &
     499        11464 :                       1.0_dp*bxy*2*drji*sinalpha*drjisinalpha)*Vgaussian*(-1.0_dp)
     500              : 
     501              :          ! Force due to gaussian term
     502        11464 :          f_nonbond(1:3, iparticle) = f_nonbond(1:3, iparticle) + dGauss(1:3)
     503              : 
     504         2866 :          IF (use_virial) THEN
     505            0 :             pv_nonbond(1, 1:3) = pv_nonbond(1, 1:3) + rji(1)*dGauss(1:3)
     506            0 :             pv_nonbond(2, 1:3) = pv_nonbond(2, 1:3) + rji(2)*dGauss(1:3)
     507            0 :             pv_nonbond(3, 1:3) = pv_nonbond(3, 1:3) + rji(3)*dGauss(1:3)
     508              :          END IF
     509              :          !END Gaussian%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     510              : 
     511         2866 :          AO = gal21%AO1 + gal21%AO2*gal21%gcn(jparticle)
     512         2866 :          BO = gal21%BO1 + gal21%BO2*gal21%gcn(jparticle)
     513              : 
     514              :          !Derivativ of the Tang and Toennies term
     515              :          dTT(:) = (-(AO*BO + (BO**7)*gal21%c/720)*EXP(-BO*drji) + 6*(gal21%c/drji**7)* &
     516              :                    (1.0 - EXP(-BO*drji) &
     517              :                     - BO*drji*EXP(-BO*drji) &
     518              :                     - (((BO*drji)**2)/2)*EXP(-BO*drji) &
     519              :                     - (((BO*drji)**3)/6)*EXP(-BO*drji) &
     520              :                     - (((BO*drji)**4)/24)*EXP(-BO*drji) &
     521              :                     - (((BO*drji)**5)/120)*EXP(-BO*drji) &
     522              :                     - (((BO*drji)**6)/720)*EXP(-BO*drji)) &
     523        11464 :                    )*rji_hat(:)
     524              : 
     525              :          ! Force of Tang & Toennies
     526        11464 :          f_nonbond(1:3, iparticle) = f_nonbond(1:3, iparticle) - dTT(1:3)
     527              : 
     528         2866 :          IF (use_virial) THEN
     529            0 :             pv_nonbond(1, 1:3) = pv_nonbond(1, 1:3) - rji(1)*dTT(1:3)
     530            0 :             pv_nonbond(2, 1:3) = pv_nonbond(2, 1:3) - rji(2)*dTT(1:3)
     531            0 :             pv_nonbond(3, 1:3) = pv_nonbond(3, 1:3) - rji(3)*dTT(1:3)
     532              :          END IF
     533              : 
     534              :       END IF
     535              : 
     536         5732 :    END SUBROUTINE gal21_forces
     537              : 
     538              : ! **************************************************************************************************
     539              : !> \brief Derivativ of the second component of angular dependance
     540              : !> \param gal21 ...
     541              : !> \param r_last_update_pbc ...
     542              : !> \param iparticle ...
     543              : !> \param jparticle ...
     544              : !> \param f_nonbond ...
     545              : !> \param pv_nonbond ...
     546              : !> \param use_virial ...
     547              : !> \param particle_set ...
     548              : !> \param cell ...
     549              : !> \param anglepart ...
     550              : !> \param sum_weight ...
     551              : ! **************************************************************************************************
     552         2866 :    SUBROUTINE somme_d(gal21, r_last_update_pbc, iparticle, jparticle, &
     553         2866 :                       f_nonbond, pv_nonbond, use_virial, particle_set, cell, anglepart, sum_weight)
     554              :       TYPE(gal21_pot_type), POINTER                      :: gal21
     555              :       TYPE(pos_type), DIMENSION(:), POINTER              :: r_last_update_pbc
     556              :       INTEGER, INTENT(IN)                                :: iparticle, jparticle
     557              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT)      :: f_nonbond, pv_nonbond
     558              :       LOGICAL, INTENT(IN)                                :: use_virial
     559              :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     560              :       TYPE(cell_type), POINTER                           :: cell
     561              :       REAL(KIND=dp), INTENT(IN)                          :: anglepart, sum_weight
     562              : 
     563              :       CHARACTER(LEN=2)                                   :: element_symbol_k
     564              :       INTEGER                                            :: kparticle, natom
     565              :       REAL(KIND=dp)                                      :: drki, dwdr(3), rji(3), rki(3), &
     566              :                                                             rki_hat(3), weight_rji
     567              : 
     568         2866 :       rji(:) = pbc(r_last_update_pbc(jparticle)%r(:), r_last_update_pbc(iparticle)%r(:), cell)
     569        11464 :       weight_rji = EXP(-SQRT(DOT_PRODUCT(rji, rji))/gal21%r1)
     570              : 
     571         2866 :       natom = SIZE(particle_set)
     572      2496286 :       DO kparticle = 1, natom !Loop on every atom of the system
     573              :          CALL get_atomic_kind(atomic_kind=particle_set(kparticle)%atomic_kind, &
     574      2493420 :                               element_symbol=element_symbol_k)
     575      2493420 :          IF (element_symbol_k /= gal21%met1 .AND. element_symbol_k /= gal21%met2) CYCLE !Keep only metals
     576       550272 :          rki(:) = pbc(r_last_update_pbc(kparticle)%r(:), r_last_update_pbc(iparticle)%r(:), cell)
     577      2201088 :          IF (SQRT(DOT_PRODUCT(rki, rki)) > gal21%rcutsq) CYCLE !Keep only those within cutoff distance of the oxygen atom of the evaluated pair (the omega ensemble)
     578      2201088 :          drki = SQRT(DOT_PRODUCT(rki, rki))
     579      2201088 :          rki_hat(:) = rki(:)/drki
     580              : 
     581      2201088 :          IF (element_symbol_k == gal21%met1) dwdr(:) = (-1.0_dp)*(1.0_dp/gal21%r1)*EXP(-drki/gal21%r1)*rki_hat(:) !Build the sum of derivativs
     582       550272 :          IF (element_symbol_k == gal21%met2) dwdr(:) = (-1.0_dp)*(1.0_dp/gal21%r2)*EXP(-drki/gal21%r2)*rki_hat(:)
     583              : 
     584              :          f_nonbond(1:3, iparticle) = f_nonbond(1:3, iparticle) + dwdr(1:3)*weight_rji &
     585      2201088 :                                      *weight_rji*anglepart/(sum_weight**2)
     586              : 
     587       553138 :          IF (use_virial) THEN
     588              :             pv_nonbond(1, 1:3) = pv_nonbond(1, 1:3) + rki(1)*dwdr(1:3)*weight_rji &
     589            0 :                                  *weight_rji*anglepart/(sum_weight**2)
     590              :             pv_nonbond(2, 1:3) = pv_nonbond(2, 1:3) + rki(2)*dwdr(1:3)*weight_rji &
     591            0 :                                  *weight_rji*anglepart/(sum_weight**2)
     592              :             pv_nonbond(3, 1:3) = pv_nonbond(3, 1:3) + rki(3)*dwdr(1:3)*weight_rji &
     593            0 :                                  *weight_rji*anglepart/(sum_weight**2)
     594              :          END IF
     595              : 
     596              :       END DO
     597              : 
     598         2866 :    END SUBROUTINE somme_d
     599              : 
     600              : ! **************************************************************************************************
     601              : !> \brief Derivativ of the third component of angular term
     602              : !> \param gal21 ...
     603              : !> \param r_last_update_pbc ...
     604              : !> \param iparticle ...
     605              : !> \param jparticle ...
     606              : !> \param f_nonbond ...
     607              : !> \param pv_nonbond ...
     608              : !> \param use_virial ...
     609              : !> \param prefactor ...
     610              : !> \param cell ...
     611              : !> \param particle_set ...
     612              : !> \param nvec ...
     613              : ! **************************************************************************************************
     614         2866 :    SUBROUTINE angular_d(gal21, r_last_update_pbc, iparticle, jparticle, f_nonbond, &
     615         2866 :                         pv_nonbond, use_virial, prefactor, cell, particle_set, nvec)
     616              :       TYPE(gal21_pot_type), POINTER                      :: gal21
     617              :       TYPE(pos_type), DIMENSION(:), POINTER              :: r_last_update_pbc
     618              :       INTEGER, INTENT(IN)                                :: iparticle, jparticle
     619              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT)      :: f_nonbond, pv_nonbond
     620              :       LOGICAL, INTENT(IN)                                :: use_virial
     621              :       REAL(KIND=dp), INTENT(IN)                          :: prefactor
     622              :       TYPE(cell_type), POINTER                           :: cell
     623              :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     624              :       REAL(KIND=dp), DIMENSION(3)                        :: nvec
     625              : 
     626              :       CHARACTER(LEN=2)                                   :: element_symbol
     627              :       INTEGER                                            :: count_h, iatom, index_h1, index_h2, natom
     628              :       REAL(KIND=dp)                                      :: a1, a2, a3, a4, BH, costheta, &
     629              :                                                             dsumdtheta, h_max_dist, theta
     630              :       REAL(KIND=dp), DIMENSION(3)                        :: dangular, dcostheta, rih, rih1, rih2, &
     631              :                                                             rix, rix_hat, rjh1, rjh2, rji, rji_hat
     632              : 
     633         2866 :       count_h = 0
     634         2866 :       index_h1 = 0
     635         2866 :       index_h2 = 0
     636         2866 :       h_max_dist = 2.1_dp ! 1.1 angstrom
     637         2866 :       natom = SIZE(particle_set)
     638              : 
     639      2496286 :       DO iatom = 1, natom !Loop on every atom of the system
     640              :          CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
     641      2493420 :                               element_symbol=element_symbol)
     642      2493420 :          IF (element_symbol /= "H") CYCLE !Kepp only hydrogen
     643      1295432 :          rih(:) = pbc(r_last_update_pbc(iparticle)%r(:), r_last_update_pbc(iatom)%r(:), cell)
     644      5181728 :          IF (SQRT(DOT_PRODUCT(rih, rih)) >= h_max_dist) CYCLE !Keep only hydrogen that are bounded to the considered O
     645         5732 :          count_h = count_h + 1
     646         8598 :          IF (count_h == 1) THEN
     647              :             index_h1 = iatom
     648         2866 :          ELSEIF (count_h == 2) THEN
     649         2866 :             index_h2 = iatom
     650              :          END IF
     651              :       END DO
     652              : 
     653              :       ! Abort if the oxygen is not part of a water molecule (2 H)
     654         2866 :       IF (count_h /= 2) THEN
     655              :          CALL cp_abort(__LOCATION__, &
     656            0 :                        " Error: Found "//cp_to_string(count_h)//" H atoms for O atom "//cp_to_string(iparticle))
     657              :       END IF
     658              : 
     659         2866 :       a1 = gal21%a11 + gal21%a12*gal21%gcn(jparticle) + gal21%a13*gal21%gcn(jparticle)*gal21%gcn(jparticle)
     660         2866 :       a2 = gal21%a21 + gal21%a22*gal21%gcn(jparticle) + gal21%a23*gal21%gcn(jparticle)*gal21%gcn(jparticle)
     661         2866 :       a3 = gal21%a31 + gal21%a32*gal21%gcn(jparticle) + gal21%a33*gal21%gcn(jparticle)*gal21%gcn(jparticle)
     662         2866 :       a4 = gal21%a41 + gal21%a42*gal21%gcn(jparticle) + gal21%a43*gal21%gcn(jparticle)*gal21%gcn(jparticle)
     663              : 
     664         2866 :       rji(:) = pbc(r_last_update_pbc(jparticle)%r(:), r_last_update_pbc(iparticle)%r(:), cell)
     665        20062 :       rji_hat(:) = rji(:)/SQRT(DOT_PRODUCT(rji, rji)) ! hat = pure directional component of a given vector
     666              : 
     667              :       !dipole vector rix of the H2O molecule
     668         2866 :       rih1(:) = pbc(r_last_update_pbc(iparticle)%r(:), r_last_update_pbc(index_h1)%r(:), cell)
     669         2866 :       rih2(:) = pbc(r_last_update_pbc(iparticle)%r(:), r_last_update_pbc(index_h2)%r(:), cell)
     670        11464 :       rix(:) = rih1(:) + rih2(:) ! build the dipole vector rix of the H2O molecule
     671        20062 :       rix_hat(:) = rix(:)/SQRT(DOT_PRODUCT(rix, rix)) ! hat = pure directional component of a given vector
     672        20062 :       costheta = DOT_PRODUCT(rix, nvec)/SQRT(DOT_PRODUCT(rix, rix)) ! Theta is the angle between the normal to the surface and the dipole
     673         2866 :       IF (costheta < -1.0_dp) costheta = -1.0_dp
     674         2866 :       IF (costheta > +1.0_dp) costheta = +1.0_dp
     675         2866 :       theta = ACOS(costheta) ! Theta is the angle between the normal to the surface and the dipole
     676              : 
     677              :       ! Calculation of partial derivativ of the angular components
     678              :       dsumdtheta = -1.0_dp*a1*SIN(theta) - a2*2.0_dp*SIN(2.0_dp*theta) - &
     679         2866 :                    a3*3.0_dp*SIN(3.0_dp*theta) - a4*4.0_dp*SIN(4.0_dp*theta)
     680        20062 :       dcostheta(:) = (1.0_dp/SQRT(DOT_PRODUCT(rix, rix)))*(nvec(:) - costheta*rix_hat(:))
     681        11464 :       dangular(:) = prefactor*dsumdtheta*(-1.0_dp/SIN(theta))*dcostheta(:)
     682              : 
     683              :       !Force due to the third component of the derivativ of the angular term
     684        11464 :       f_nonbond(1:3, iparticle) = f_nonbond(1:3, iparticle) - dangular(1:3)*2.0_dp !(one per H)
     685        11464 :       f_nonbond(1:3, index_h1) = f_nonbond(1:3, index_h1) + dangular(1:3)
     686        11464 :       f_nonbond(1:3, index_h2) = f_nonbond(1:3, index_h2) + dangular(1:3)
     687              : 
     688         2866 :       IF (use_virial) THEN
     689            0 :          pv_nonbond(1, 1:3) = pv_nonbond(1, 1:3) + rix(1)*dangular(1:3)
     690            0 :          pv_nonbond(2, 1:3) = pv_nonbond(2, 1:3) + rix(2)*dangular(1:3)
     691            0 :          pv_nonbond(3, 1:3) = pv_nonbond(3, 1:3) + rix(3)*dangular(1:3)
     692              :       END IF
     693              : 
     694         2866 :       BH = gal21%BH1 + gal21%gcn(jparticle)*gal21%BH2
     695              : 
     696         2866 :       rjh1(:) = pbc(r_last_update_pbc(jparticle)%r(:), r_last_update_pbc(index_h1)%r(:), cell)
     697              :       f_nonbond(1:3, index_h1) = f_nonbond(1:3, index_h1) + (gal21%AH2*gal21%gcn(jparticle) + gal21%AH1)* &
     698        20062 :                                  BH*EXP(-BH*SQRT(DOT_PRODUCT(rjh1, rjh1)))*rjh1(:)/SQRT(DOT_PRODUCT(rjh1, rjh1))
     699              : 
     700         2866 :       IF (use_virial) THEN
     701              :          pv_nonbond(1, 1:3) = pv_nonbond(1, 1:3) + rjh1(1)*((gal21%AH2*gal21%gcn(jparticle) + gal21%AH1)* &
     702              :                                                             BH*EXP(-BH*SQRT(DOT_PRODUCT(rjh1, rjh1)))) &
     703            0 :                               *rjh1(:)/SQRT(DOT_PRODUCT(rjh1, rjh1))
     704              :          pv_nonbond(2, 1:3) = pv_nonbond(2, 1:3) + rjh1(2)*((gal21%AH2*gal21%gcn(jparticle) + gal21%AH1)* &
     705              :                                                             BH*EXP(-BH*SQRT(DOT_PRODUCT(rjh1, rjh1)))) &
     706            0 :                               *rjh1(:)/SQRT(DOT_PRODUCT(rjh1, rjh1))
     707              :          pv_nonbond(3, 1:3) = pv_nonbond(3, 1:3) + rjh1(3)*((gal21%AH2*gal21%gcn(jparticle) + gal21%AH1)* &
     708              :                                                             BH*EXP(-BH*SQRT(DOT_PRODUCT(rjh1, rjh1)))) &
     709            0 :                               *rjh1(:)/SQRT(DOT_PRODUCT(rjh1, rjh1))
     710              :       END IF
     711              : 
     712         2866 :       rjh2(:) = pbc(r_last_update_pbc(jparticle)%r(:), r_last_update_pbc(index_h2)%r(:), cell)
     713              :       f_nonbond(1:3, index_h2) = f_nonbond(1:3, index_h2) + ((gal21%AH2*gal21%gcn(jparticle) + gal21%AH1)* &
     714              :                                                              BH*EXP(-BH*SQRT(DOT_PRODUCT(rjh2, rjh2)))) &
     715        20062 :                                  *rjh2(:)/SQRT(DOT_PRODUCT(rjh2, rjh2))
     716              : 
     717         2866 :       IF (use_virial) THEN
     718              :          pv_nonbond(1, 1:3) = pv_nonbond(1, 1:3) + rjh2(1)*((gal21%AH2*gal21%gcn(jparticle) + gal21%AH1)* &
     719              :                                                             BH*EXP(-BH*SQRT(DOT_PRODUCT(rjh2, rjh2)))) &
     720            0 :                               *rjh2(:)/SQRT(DOT_PRODUCT(rjh2, rjh2))
     721              :          pv_nonbond(2, 1:3) = pv_nonbond(2, 1:3) + rjh2(2)*((gal21%AH2*gal21%gcn(jparticle) + gal21%AH1)* &
     722              :                                                             BH*EXP(-BH*SQRT(DOT_PRODUCT(rjh2, rjh2)))) &
     723            0 :                               *rjh2(:)/SQRT(DOT_PRODUCT(rjh2, rjh2))
     724              :          pv_nonbond(3, 1:3) = pv_nonbond(3, 1:3) + rjh2(3)*((gal21%AH2*gal21%gcn(jparticle) + gal21%AH1)* &
     725              :                                                             BH*EXP(-BH*SQRT(DOT_PRODUCT(rjh2, rjh2)))) &
     726            0 :                               *rjh2(:)/SQRT(DOT_PRODUCT(rjh2, rjh2))
     727              :       END IF
     728              : 
     729         2866 :    END SUBROUTINE angular_d
     730              : 
     731              : ! **************************************************************************************************
     732              : !> \brief ...
     733              : !> \param nonbonded ...
     734              : !> \param potparm ...
     735              : !> \param glob_loc_list ...
     736              : !> \param glob_cell_v ...
     737              : !> \param glob_loc_list_a ...
     738              : !> \param cell ...
     739              : !> \par History
     740              : ! **************************************************************************************************
     741            2 :    SUBROUTINE setup_gal21_arrays(nonbonded, potparm, glob_loc_list, glob_cell_v, &
     742              :                                  glob_loc_list_a, cell)
     743              :       TYPE(fist_neighbor_type), POINTER                  :: nonbonded
     744              :       TYPE(pair_potential_pp_type), POINTER              :: potparm
     745              :       INTEGER, DIMENSION(:, :), POINTER                  :: glob_loc_list
     746              :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: glob_cell_v
     747              :       INTEGER, DIMENSION(:), POINTER                     :: glob_loc_list_a
     748              :       TYPE(cell_type), POINTER                           :: cell
     749              : 
     750              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'setup_gal21_arrays'
     751              : 
     752              :       INTEGER                                            :: handle, i, iend, igrp, ikind, ilist, &
     753              :                                                             ipair, istart, jkind, nkinds, npairs, &
     754              :                                                             npairs_tot
     755            2 :       INTEGER, DIMENSION(:), POINTER                     :: work_list, work_list2
     756            2 :       INTEGER, DIMENSION(:, :), POINTER                  :: list
     757              :       REAL(KIND=dp), DIMENSION(3)                        :: cell_v, cvi
     758            2 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: rwork_list
     759              :       TYPE(neighbor_kind_pairs_type), POINTER            :: neighbor_kind_pair
     760              :       TYPE(pair_potential_single_type), POINTER          :: pot
     761              : 
     762            0 :       CPASSERT(.NOT. ASSOCIATED(glob_loc_list))
     763            2 :       CPASSERT(.NOT. ASSOCIATED(glob_loc_list_a))
     764            2 :       CPASSERT(.NOT. ASSOCIATED(glob_cell_v))
     765            2 :       CALL timeset(routineN, handle)
     766            2 :       npairs_tot = 0
     767            2 :       nkinds = SIZE(potparm%pot, 1)
     768           56 :       DO ilist = 1, nonbonded%nlists
     769           54 :          neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist)
     770           54 :          npairs = neighbor_kind_pair%npairs
     771           54 :          IF (npairs == 0) CYCLE
     772          336 :          Kind_Group_Loop1: DO igrp = 1, neighbor_kind_pair%ngrp_kind
     773          316 :             istart = neighbor_kind_pair%grp_kind_start(igrp)
     774          316 :             iend = neighbor_kind_pair%grp_kind_end(igrp)
     775          316 :             ikind = neighbor_kind_pair%ij_kind(1, igrp)
     776          316 :             jkind = neighbor_kind_pair%ij_kind(2, igrp)
     777          316 :             pot => potparm%pot(ikind, jkind)%pot
     778          316 :             npairs = iend - istart + 1
     779          316 :             IF (pot%no_mb) CYCLE
     780           90 :             DO i = 1, SIZE(pot%type)
     781          334 :                IF (pot%type(i) == gal21_type) npairs_tot = npairs_tot + npairs
     782              :             END DO
     783              :          END DO Kind_Group_Loop1
     784              :       END DO
     785            6 :       ALLOCATE (work_list(npairs_tot))
     786            4 :       ALLOCATE (work_list2(npairs_tot))
     787            6 :       ALLOCATE (glob_loc_list(2, npairs_tot))
     788            6 :       ALLOCATE (glob_cell_v(3, npairs_tot))
     789              :       ! Fill arrays with data
     790            2 :       npairs_tot = 0
     791           56 :       DO ilist = 1, nonbonded%nlists
     792           54 :          neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist)
     793           54 :          npairs = neighbor_kind_pair%npairs
     794           54 :          IF (npairs == 0) CYCLE
     795          336 :          Kind_Group_Loop2: DO igrp = 1, neighbor_kind_pair%ngrp_kind
     796          316 :             istart = neighbor_kind_pair%grp_kind_start(igrp)
     797          316 :             iend = neighbor_kind_pair%grp_kind_end(igrp)
     798          316 :             ikind = neighbor_kind_pair%ij_kind(1, igrp)
     799          316 :             jkind = neighbor_kind_pair%ij_kind(2, igrp)
     800          316 :             list => neighbor_kind_pair%list
     801         1264 :             cvi = neighbor_kind_pair%cell_vector
     802          316 :             pot => potparm%pot(ikind, jkind)%pot
     803          316 :             npairs = iend - istart + 1
     804          316 :             IF (pot%no_mb) CYCLE
     805          234 :             cell_v = MATMUL(cell%hmat, cvi)
     806           90 :             DO i = 1, SIZE(pot%type)
     807              :                ! gal21
     808          334 :                IF (pot%type(i) == gal21_type) THEN
     809        17618 :                   DO ipair = 1, npairs
     810       105600 :                      glob_loc_list(:, npairs_tot + ipair) = list(:, istart - 1 + ipair)
     811        70418 :                      glob_cell_v(1:3, npairs_tot + ipair) = cell_v(1:3)
     812              :                   END DO
     813           18 :                   npairs_tot = npairs_tot + npairs
     814              :                END IF
     815              :             END DO
     816              :          END DO Kind_Group_Loop2
     817              :       END DO
     818              :       ! Order the arrays w.r.t. the first index of glob_loc_list
     819            2 :       CALL sort(glob_loc_list(1, :), npairs_tot, work_list)
     820        17602 :       DO ipair = 1, npairs_tot
     821        17602 :          work_list2(ipair) = glob_loc_list(2, work_list(ipair))
     822              :       END DO
     823        35204 :       glob_loc_list(2, :) = work_list2
     824            2 :       DEALLOCATE (work_list2)
     825            6 :       ALLOCATE (rwork_list(3, npairs_tot))
     826        17602 :       DO ipair = 1, npairs_tot
     827       140802 :          rwork_list(:, ipair) = glob_cell_v(:, work_list(ipair))
     828              :       END DO
     829       140804 :       glob_cell_v = rwork_list
     830            2 :       DEALLOCATE (rwork_list)
     831            2 :       DEALLOCATE (work_list)
     832            6 :       ALLOCATE (glob_loc_list_a(npairs_tot))
     833        35204 :       glob_loc_list_a = glob_loc_list(1, :)
     834            2 :       CALL timestop(handle)
     835            4 :    END SUBROUTINE setup_gal21_arrays
     836              : 
     837              : ! **************************************************************************************************
     838              : !> \brief ...
     839              : !> \param glob_loc_list ...
     840              : !> \param glob_cell_v ...
     841              : !> \param glob_loc_list_a ...
     842              : ! **************************************************************************************************
     843            1 :    SUBROUTINE destroy_gal21_arrays(glob_loc_list, glob_cell_v, glob_loc_list_a)
     844              :       INTEGER, DIMENSION(:, :), POINTER                  :: glob_loc_list
     845              :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: glob_cell_v
     846              :       INTEGER, DIMENSION(:), POINTER                     :: glob_loc_list_a
     847              : 
     848            1 :       IF (ASSOCIATED(glob_loc_list)) THEN
     849            1 :          DEALLOCATE (glob_loc_list)
     850              :       END IF
     851            1 :       IF (ASSOCIATED(glob_loc_list_a)) THEN
     852            1 :          DEALLOCATE (glob_loc_list_a)
     853              :       END IF
     854            1 :       IF (ASSOCIATED(glob_cell_v)) THEN
     855            1 :          DEALLOCATE (glob_cell_v)
     856              :       END IF
     857              : 
     858            1 :    END SUBROUTINE destroy_gal21_arrays
     859              : 
     860              : ! **************************************************************************************************
     861              : !> \brief prints the number of OH- ions or H3O+ ions near surface
     862              : !> \param nr_ions number of ions
     863              : !> \param mm_section ...
     864              : !> \param para_env ...
     865              : !> \param print_oh flag indicating if number OH- is printed
     866              : !> \param print_h3o flag indicating if number H3O+ is printed
     867              : !> \param print_o flag indicating if number O^(2-) is printed
     868              : ! **************************************************************************************************
     869            0 :    SUBROUTINE print_nr_ions_gal21(nr_ions, mm_section, para_env, print_oh, &
     870              :                                   print_h3o, print_o)
     871              :       INTEGER, INTENT(INOUT)                             :: nr_ions
     872              :       TYPE(section_vals_type), POINTER                   :: mm_section
     873              :       TYPE(mp_para_env_type), OPTIONAL, POINTER          :: para_env
     874              :       LOGICAL, INTENT(IN)                                :: print_oh, print_h3o, print_o
     875              : 
     876              :       INTEGER                                            :: iw
     877              :       TYPE(cp_logger_type), POINTER                      :: logger
     878              : 
     879            0 :       NULLIFY (logger)
     880              : 
     881            0 :       CALL para_env%sum(nr_ions)
     882            0 :       logger => cp_get_default_logger()
     883              : 
     884              :       iw = cp_print_key_unit_nr(logger, mm_section, "PRINT%PROGRAM_RUN_INFO", &
     885            0 :                                 extension=".mmLog")
     886              : 
     887            0 :       IF (iw > 0 .AND. nr_ions > 0 .AND. print_oh) THEN
     888            0 :          WRITE (iw, '(/,A,T71,I10,/)') " gal21: number of OH- ions at surface", nr_ions
     889              :       END IF
     890            0 :       IF (iw > 0 .AND. nr_ions > 0 .AND. print_h3o) THEN
     891            0 :          WRITE (iw, '(/,A,T71,I10,/)') " gal21: number of H3O+ ions at surface", nr_ions
     892              :       END IF
     893            0 :       IF (iw > 0 .AND. nr_ions > 0 .AND. print_o) THEN
     894            0 :          WRITE (iw, '(/,A,T71,I10,/)') " gal21: number of O^2- ions at surface", nr_ions
     895              :       END IF
     896              : 
     897            0 :       CALL cp_print_key_finished_output(iw, logger, mm_section, "PRINT%PROGRAM_RUN_INFO")
     898              : 
     899            0 :    END SUBROUTINE print_nr_ions_gal21
     900              : 
     901              : END MODULE manybody_gal21
        

Generated by: LCOV version 2.0-1