Skip to content

Commit

Permalink
Interface to GRRM17 code (#715)
Browse files Browse the repository at this point in the history 
* GRRM interface

* Energy correction for NDDO methods

* Protect from Debug runs, adjust regtest value

* prettify/doxify run

Co-authored-by: Tiziano Müller <tm@dev-zero.ch>
  • Loading branch information
@juerghutter @dev-zero
juerghutter and dev-zero committed Jan 8, 2020
1 parent 502ab3a commit cf243b8
Show file tree
Hide file tree
Showing 7 changed files with 185 additions and 13 deletions.
26 changes: 22 additions & 4 deletions src/force_env_methods.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -111,6 +111,9 @@ MODULE force_env_methods
mixed_environment_type
USE mixed_environment_utils, ONLY: get_subsys_map_index,&
mixed_map_forces
USE molecule_kind_list_types, ONLY: molecule_kind_list_type
USE molecule_kind_types, ONLY: get_molecule_kind,&
molecule_kind_type
USE optimize_dmfet_potential, ONLY: build_full_dm,&
check_dmfet,&
prepare_dmfet_opt,&
Expand Down Expand Up @@ -208,8 +211,9 @@ RECURSIVE SUBROUTINE force_env_calc_energy_force(force_env, calc_force, &
routineP = moduleN//':'//routineN
REAL(kind=dp), PARAMETER :: ateps = 1.0E-6_dp

INTEGER :: i, j, nat, ndigits, output_unit, &
print_forces, print_grrm
INTEGER :: i, ikind, j, nat, ndigits, nfixed_atoms, &
nfixed_atoms_total, nkind, &
output_unit, print_forces, print_grrm
LOGICAL :: calculate_forces, &
calculate_stress_tensor, &
energy_consistency, eval_ef, &
Expand All @@ -222,6 +226,9 @@ RECURSIVE SUBROUTINE force_env_calc_energy_force(force_env, calc_force, &
TYPE(cell_type), POINTER :: cell
TYPE(cp_logger_type), POINTER :: logger
TYPE(cp_subsys_type), POINTER :: subsys
TYPE(molecule_kind_list_type), POINTER :: molecule_kinds
TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
TYPE(molecule_kind_type), POINTER :: molecule_kind
TYPE(particle_list_type), POINTER :: core_particles, particles, &
shell_particles
TYPE(virial_type), POINTER :: virial
Expand Down Expand Up @@ -435,8 +442,19 @@ RECURSIVE SUBROUTINE force_env_calc_energy_force(force_env, calc_force, &
file_position="REWIND", extension=".rrm")
IF (print_grrm > 0) THEN
CALL force_env_get(force_env, subsys=subsys)
CALL cp_subsys_get(subsys, particles=particles)
CALL write_grrm(print_grrm, particles%els, e_pot)
CALL cp_subsys_get(subsys=subsys, particles=particles, &
molecule_kinds=molecule_kinds)
! Count the number of fixed atoms
nfixed_atoms_total = 0
nkind = molecule_kinds%n_els
molecule_kind_set => molecule_kinds%els
DO ikind = 1, nkind
molecule_kind => molecule_kind_set(ikind)
CALL get_molecule_kind(molecule_kind, nfixd=nfixed_atoms)
nfixed_atoms_total = nfixed_atoms_total + nfixed_atoms
END DO
!
CALL write_grrm(print_grrm, force_env, particles%els, e_pot, fixed_atoms=nfixed_atoms_total)
END IF
CALL cp_print_key_finished_output(print_grrm, logger, force_env%force_env_section, "PRINT%GRRM")
Expand Down
34 changes: 29 additions & 5 deletions src/grrm_utils.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,9 +9,13 @@
! **************************************************************************************************
MODULE grrm_utils

USE cp_control_types, ONLY: dft_control_type
USE force_env_types, ONLY: force_env_type
USE kinds, ONLY: dp
USE particle_types, ONLY: particle_type
USE physcon, ONLY: angstrom
USE qs_energy_types, ONLY: qs_energy_type
USE qs_environment_types, ONLY: get_qs_env
#include "./base/base_uses.f90"

IMPLICIT NONE
Expand All @@ -30,44 +34,65 @@ MODULE grrm_utils
!> \brief Write GRRM interface file
!>
!> \param iounit ...
!> \param force_env ...
!> \param particles ...
!> \param energy ...
!> \param dipole ...
!> \param hessian ...
!> \param dipder ...
!> \param polar ...
!> \param fixed_atoms ...
! **************************************************************************************************
SUBROUTINE write_grrm(iounit, particles, energy, dipole, hessian, dipder, polar)
SUBROUTINE write_grrm(iounit, force_env, particles, energy, dipole, hessian, dipder, polar, &
fixed_atoms)

INTEGER, INTENT(IN) :: iounit
TYPE(force_env_type), POINTER :: force_env
TYPE(particle_type), DIMENSION(:), INTENT(IN) :: particles
REAL(KIND=dp), INTENT(IN) :: energy
REAL(KIND=dp), DIMENSION(3), INTENT(IN), OPTIONAL :: dipole
REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
OPTIONAL :: hessian, dipder
REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN), &
OPTIONAL :: polar
INTEGER, INTENT(IN), OPTIONAL :: fixed_atoms

CHARACTER(LEN=*), PARAMETER :: routineN = 'write_grrm', routineP = moduleN//':'//routineN
REAL(KIND=dp), PARAMETER :: zero = 0.0_dp

INTEGER :: i, j, natom, nc
LOGICAL :: nddo
REAL(KIND=dp) :: eout
REAL(KIND=dp), DIMENSION(5) :: fz
TYPE(dft_control_type), POINTER :: dft_control
TYPE(qs_energy_type), POINTER :: qs_energy

IF (iounit > 0) THEN
! the units depend on the qs method!
CPASSERT(ASSOCIATED(force_env%qs_env))
CALL get_qs_env(force_env%qs_env, dft_control=dft_control)
nddo = dft_control%qs_control%semi_empirical
IF (nddo) THEN
CALL get_qs_env(force_env%qs_env, energy=qs_energy)
eout = energy + qs_energy%core_self
ELSE
eout = energy
END IF
!
natom = SIZE(particles)
IF (PRESENT(fixed_atoms)) natom = natom - fixed_atoms
WRITE (iounit, "(A7)") "RESULTS"
WRITE (iounit, "(A18)") "CURRENT COORDINATE"
DO i = 1, natom
WRITE (iounit, "(A,3F24.12)") TRIM(ADJUSTL(particles(i)%atomic_kind%element_symbol)), &
particles(i)%r(1:3)*angstrom
END DO
WRITE (iounit, "(A8,3F18.12)") "ENERGY =", energy, zero, zero
WRITE (iounit, "(A8,3F18.12)") "ENERGY =", eout, zero, zero
WRITE (iounit, "(A8,3F18.12)") " =", zero, zero, zero
WRITE (iounit, "(A8,F18.12)") "S**2 =", zero
WRITE (iounit, "(A8)") "GRADIENT"
DO i = 1, natom
WRITE (iounit, "(F17.12)") particles(i)%f(1:3)
WRITE (iounit, "(F17.12)") - particles(i)%f(1:3)
END DO
IF (PRESENT(dipole)) THEN
WRITE (iounit, "(A8,3F18.12)") "DIPOLE =", dipole(1:3)
Expand All @@ -77,8 +102,7 @@ SUBROUTINE write_grrm(iounit, particles, energy, dipole, hessian, dipder, polar)
fz = zero
WRITE (iounit, "(A7)") "HESSIAN"
IF (PRESENT(hessian)) THEN
nc = SIZE(hessian, 1)
CPASSERT(nc == 3*natom)
nc = 3*natom
DO i = 1, nc, 5
DO j = i, nc
WRITE (iounit, "(5(F13.9,1X))") hessian(j, i:MIN(j, i + 4))
Expand Down
2 changes: 1 addition & 1 deletion src/input_cp2k_force_eval.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -370,7 +370,7 @@ SUBROUTINE create_f_env_print_section(section)

CALL cp_print_key_section_create(print_key, __LOCATION__, "GRRM", &
description="Controls the printing of the GRRM interface file", &
print_level=debug_print_level, filename="CP2K_GRMM")
print_level=debug_print_level + 1, filename="CP2K_GRMM")
CALL section_add_subsection(section, print_key)
CALL section_release(print_key)

Expand Down
8 changes: 5 additions & 3 deletions src/motion/vibrational_analysis.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -105,8 +105,8 @@ SUBROUTINE vb_anal(input, input_declaration, para_env, globenv)
CHARACTER(LEN=default_string_length) :: description
INTEGER :: handle, i, icoord, icoordm, icoordp, ierr, imap, ip1, ip2, iparticle1, &
iparticle2, iseq, iw, j, k, natoms, ncoord, nrep, nres, nRotTrM, nvib, output_unit, &
output_unit_eig, prep, print_grrm, proc_dist_type
iparticle2, iseq, iw, j, k, natoms, ncoord, nfrozen, nrep, nres, nRotTrM, nvib, &
output_unit, output_unit_eig, prep, print_grrm, proc_dist_type
INTEGER, DIMENSION(:), POINTER :: Clist, Mlist
LOGICAL :: calc_intens, calc_thchdata, &
do_mode_tracking, keep_rotations, &
Expand Down Expand Up @@ -414,7 +414,9 @@ SUBROUTINE vb_anal(input, input_declaration, para_env, globenv)
Hint1(j, i) = Hessian(j, i)*rmass(i)*rmass(j)*1.0E-6_dp
END DO
END DO
CALL write_grrm(print_grrm, particles, minimum_energy, hessian=Hint1)
nfrozen = SIZE(particles) - natoms
CALL write_grrm(print_grrm, f_env%force_env, particles, minimum_energy, &
hessian=Hint1, fixed_atoms=nfrozen)
DEALLOCATE (Hint1, rmass)
END IF
CALL cp_print_key_finished_output(print_grrm, logger, force_env_section, "PRINT%GRRM")
Expand Down
60 changes: 60 additions & 0 deletions tests/QS/regtest-gpw-9/C2H4_frozen_GRRM.inp
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,60 @@
&FORCE_EVAL
METHOD Quickstep
&PRINT
&GRRM
FILENAME =grrm.rrm
&END
&END
&DFT
BASIS_SET_FILE_NAME BASIS_SET
POTENTIAL_FILE_NAME POTENTIAL
&MGRID
CUTOFF 50
&END MGRID
&QS
&END QS
&SCF
MAX_SCF 5
&OT ON
&END OT
&END SCF
&XC
&XC_FUNCTIONAL PADE
&END XC_FUNCTIONAL
&END XC
&END DFT
&SUBSYS
&CELL
ABC 6.0 6.0 8.0
&END CELL
&COORD
C -0.2558206925 -0.0449868417 0.6850619693
C -0.2214915625 -0.0386665897 -0.6514432536
H -0.4100752923 0.8714166677 1.2581094627
H -0.1401850725 -0.9650616559 1.2579495323
H -0.0795226532 -0.9560364450 -1.2252531634
H -0.3538336094 0.8817708920 -1.2228724740
&END COORD
&KIND H
BASIS_SET SZV-GTH-PADE
POTENTIAL GTH-PADE-q1
&END KIND
&KIND C
BASIS_SET SZV-GTH-PADE
POTENTIAL GTH-PADE-q4
&END KIND
&END SUBSYS
&END FORCE_EVAL
&GLOBAL
PROJECT C2H4
RUN_TYPE ENERGY_FORCE
PRINT_LEVEL LOW
&END GLOBAL
&MOTION
&CONSTRAINT
&FIXED_ATOMS
LIST 3..6
&END
&END
&END

2 changes: 2 additions & 0 deletions tests/QS/regtest-gpw-9/TEST_FILES
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,4 +2,6 @@
nlcc-mix.inp 1 1.0E-12 -18.109625525665670
C2H4_GRRM.inp 1 1.0E-12 -13.924538806910600
h2o_GRRM.inp 8 8.0E-06 1465.672
C2H4_frozen_GRRM.inp 1 1.0E-12 -13.924538806910600
h2o_frozen_GRRM.inp 8 8.0E-06 872.278
#EOF
66 changes: 66 additions & 0 deletions tests/QS/regtest-gpw-9/h2o_frozen_GRRM.inp
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,66 @@
&FORCE_EVAL
METHOD Quickstep
&PRINT
&GRRM
FILENAME =grrm.rrm
&END
&END
&DFT
&QS
METHOD xTB
&END QS
&SCF
SCF_GUESS MOPAC
&OT
PRECONDITIONER FULL_SINGLE_INVERSE
MINIMIZER DIIS
&END
&OUTER_SCF
MAX_SCF 10
EPS_SCF 1.0E-8
&END
MAX_SCF 10
EPS_SCF 1.0E-8
&END SCF
&PRINT
&MOMENTS ON
PERIODIC .FALSE.
REFERENCE COM
&END
&END
&END DFT
&SUBSYS
&CELL
ABC [angstrom] 6 6 6
PERIODIC NONE
&END
&COORD
O 0.000000 0.000000 -0.065587
H 0.000000 -0.757136 0.520545
H 0.000000 0.757136 0.520545
&END COORD
&TOPOLOGY
&CENTER_COORDINATES
&END
&END
&END SUBSYS
&END FORCE_EVAL
&GLOBAL
PRINT_LEVEL LOW
PROJECT h2oVib
RUN_TYPE VIBRATIONAL_ANALYSIS
&END GLOBAL

&VIBRATIONAL_ANALYSIS
INTENSITIES .TRUE.
DX 0.001
&END

&MOTION
&CONSTRAINT
&FIXED_ATOMS
LIST 3
&END
&END
&END

0 comments on commit cf243b8

Please sign in to comment.