Skip to content
Permalink
Browse files

Scine interface (#730)

* SCINE interface (minimal implementation)

* Adjust regtest value
  • Loading branch information
juerghutter committed Jan 13, 2020
1 parent 922f2e1 commit 688071a2fdf2595b47829df0f5649b5af9d3738f
@@ -159,6 +159,7 @@ MODULE force_env_methods
USE qs_rho_types, ONLY: qs_rho_get,&
qs_rho_type
USE restraint, ONLY: restraint_control
USE scine_utils, ONLY: write_scine
USE string_utilities, ONLY: compress
USE virial_methods, ONLY: write_stress_components
USE virial_types, ONLY: cp_virial,&
@@ -213,7 +214,8 @@ RECURSIVE SUBROUTINE force_env_calc_energy_force(force_env, calc_force, &

INTEGER :: i, ikind, j, nat, ndigits, nfixed_atoms, &
nfixed_atoms_total, nkind, &
output_unit, print_forces, print_grrm
output_unit, print_forces, print_grrm, &
print_scine
LOGICAL :: calculate_forces, &
calculate_stress_tensor, &
energy_consistency, eval_ef, &
@@ -458,6 +460,16 @@ RECURSIVE SUBROUTINE force_env_calc_energy_force(force_env, calc_force, &
END IF
CALL cp_print_key_finished_output(print_grrm, logger, force_env%force_env_section, "PRINT%GRRM")
print_scine = cp_print_key_unit_nr(logger, force_env%force_env_section, "PRINT%SCINE", &
file_position="REWIND", extension=".scine")
IF (print_scine > 0) THEN
CALL force_env_get(force_env, subsys=subsys)
CALL cp_subsys_get(subsys=subsys, particles=particles)
!
CALL write_scine(print_scine, force_env, particles%els, e_pot)
END IF
CALL cp_print_key_finished_output(print_scine, logger, force_env%force_env_section, "PRINT%SCINE")
! Atomic stress
output_unit = cp_print_key_unit_nr(logger, force_env%force_env_section, "PRINT%PROGRAM_RUN_INFO", &
extension=".Log")
@@ -370,7 +370,13 @@ 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 + 1, filename="CP2K_GRMM")
print_level=debug_print_level + 1, filename="CP2K_GRRM")
CALL section_add_subsection(section, print_key)
CALL section_release(print_key)

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

@@ -74,6 +74,7 @@ MODULE vibrational_analysis
rep_env_create
USE replica_types, ONLY: rep_env_release,&
replica_env_type
USE scine_utils, ONLY: write_scine
USE util, ONLY: sort
#include "../base/base_uses.f90"
@@ -106,7 +107,7 @@ 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, nfrozen, nrep, nres, nRotTrM, nvib, &
output_unit, output_unit_eig, prep, print_grrm, proc_dist_type
output_unit, output_unit_eig, prep, print_grrm, print_scine, proc_dist_type
INTEGER, DIMENSION(:), POINTER :: Clist, Mlist
LOGICAL :: calc_intens, calc_thchdata, &
do_mode_tracking, keep_rotations, &
@@ -421,6 +422,20 @@ SUBROUTINE vb_anal(input, input_declaration, para_env, globenv)
END IF
CALL cp_print_key_finished_output(print_grrm, logger, force_env_section, "PRINT%GRRM")
!
! Print SCINE interface file
print_scine = cp_print_key_unit_nr(logger, force_env_section, "PRINT%SCINE", &
file_position="REWIND", extension=".scine")
IF (print_scine > 0) THEN
DO i = 1, natoms
imap = Mlist(i)
particles(imap)%f(1:3) = rep_env%f((imap - 1)*3 + 1:(imap - 1)*3 + 3, 1)
ENDDO
nfrozen = SIZE(particles) - natoms
CPASSERT(nfrozen == 0)
CALL write_scine(print_scine, f_env%force_env, particles, minimum_energy, hessian=Hessian)
END IF
CALL cp_print_key_finished_output(print_scine, logger, force_env_section, "PRINT%SCINE")
!
nvib = ncoord - nRotTrM
ALLOCATE (H_eigval1(ncoord))
ALLOCATE (H_eigval2(SIZE(D, 2)))
@@ -0,0 +1,171 @@
!--------------------------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 - 2020 CP2K developers group !
!--------------------------------------------------------------------------------------------------!

! **************************************************************************************************
!> \brief SCINE interface
!> \author JGH - 01.2020
! **************************************************************************************************
MODULE scine_utils

USE cell_types, ONLY: cell_type
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

PRIVATE

CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'scine_utils'

PUBLIC :: write_scine

! **************************************************************************************************

CONTAINS

! **************************************************************************************************
!> \brief Write SCINE interface file
!>
!> \param iounit ...
!> \param force_env ...
!> \param particles ...
!> \param energy ...
!> \param hessian ...
! **************************************************************************************************
SUBROUTINE write_scine(iounit, force_env, particles, energy, hessian)

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(:, :), INTENT(IN), &
OPTIONAL :: hessian

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

CHARACTER(LEN=20) :: sunit
INTEGER :: i, j, natom, nc
LOGICAL :: nddo
REAL(KIND=dp) :: eout
REAL(KIND=dp), DIMENSION(5) :: fz
TYPE(cell_type), POINTER :: cell
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
!
CALL get_qs_env(force_env%qs_env, cell=cell)
natom = SIZE(particles)
WRITE (iounit, "(A,A)") "title: ", "'SCINE Interface File'"
sunit = "'hartree'"
WRITE (iounit, "(A,F18.12,A,A)") "energy: [", eout, ", "//TRIM(sunit), " ]"
sunit = "'angstrom'"
WRITE (iounit, "(A)") "system:"
WRITE (iounit, "(T2,A,A)") "unit: ", TRIM(sunit)
WRITE (iounit, "(T2,A)") "cell: "
WRITE (iounit, "(T4,A,F24.12,',',F24.12,',',F24.12,A)") &
"- A: [", cell%hmat(1:3, 1)*angstrom, " ]"
WRITE (iounit, "(T4,A,F24.12,',',F24.12,',',F24.12,A)") &
"- B: [", cell%hmat(1:3, 2)*angstrom, " ]"
WRITE (iounit, "(T4,A,F24.12,',',F24.12,',',F24.12,A)") &
"- C: [", cell%hmat(1:3, 3)*angstrom, " ]"
WRITE (iounit, "(T2,A,L1,', ',L1,', ',L1,' ]')") "periodicity: [ ", (cell%perd == 1)
WRITE (iounit, "(T2,A)") "coordinates: "
DO i = 1, natom
WRITE (iounit, "(T4,A,A2,A,F24.12,',',F24.12,',',F24.12,A)") &
"- ", TRIM(ADJUSTL(particles(i)%atomic_kind%element_symbol))//":", &
" [", particles(i)%r(1:3)*angstrom, " ]"
END DO
WRITE (iounit, "(A)") "gradient:"
sunit = "'hartree/bohr'"
WRITE (iounit, "(T2,A,A)") "unit: ", TRIM(sunit)
WRITE (iounit, "(T2,A)") "values: "
DO i = 1, natom
WRITE (iounit, "(T4,A,A2,A,F24.12,',',F24.12,',',F24.12,A)") &
"- ", TRIM(ADJUSTL(particles(i)%atomic_kind%element_symbol))//":", &
" [", particles(i)%f(1:3), " ]"
END DO
fz = zero
IF (PRESENT(hessian)) THEN
WRITE (iounit, "(A)") "hessian:"
sunit = "'hartree/bohr^2'"
WRITE (iounit, "(T2,A,A)") "unit: ", TRIM(sunit)
DO i = 1, natom
WRITE (iounit, "(T4,A)") TRIM(ADJUSTL(particles(i)%atomic_kind%element_symbol))//":"
nc = 3*(i - 1) + 1
WRITE (iounit, "(T6,'X: [')")
DO j = 1, nc, 5
IF (nc > j + 4) THEN
WRITE (iounit, "(T12,5(F20.12,','))") hessian(nc, j:j + 4)
ELSEIF (nc == j + 4) THEN
WRITE (iounit, "(T12,4(F20.12,','),F20.12,' ]')") hessian(nc, j:nc)
ELSEIF (nc == j + 3) THEN
WRITE (iounit, "(T12,3(F20.12,','),F20.12,' ]')") hessian(nc, j:nc)
ELSEIF (nc == j + 2) THEN
WRITE (iounit, "(T12,2(F20.12,','),F20.12,' ]')") hessian(nc, j:nc)
ELSEIF (nc == j + 1) THEN
WRITE (iounit, "(T12,1(F20.12,','),F20.12,' ]')") hessian(nc, j:nc)
ELSEIF (nc == j) THEN
WRITE (iounit, "(T12,F20.12,' ]')") hessian(nc, j:nc)
END IF
END DO
nc = 3*(i - 1) + 2
WRITE (iounit, "(T6,'Y: [')")
DO j = 1, nc, 5
IF (nc > j + 4) THEN
WRITE (iounit, "(T12,5(F20.12,','))") hessian(nc, j:j + 4)
ELSEIF (nc == j + 4) THEN
WRITE (iounit, "(T12,4(F20.12,','),F20.12,' ]')") hessian(nc, j:nc)
ELSEIF (nc == j + 3) THEN
WRITE (iounit, "(T12,3(F20.12,','),F20.12,' ]')") hessian(nc, j:nc)
ELSEIF (nc == j + 2) THEN
WRITE (iounit, "(T12,2(F20.12,','),F20.12,' ]')") hessian(nc, j:nc)
ELSEIF (nc == j + 1) THEN
WRITE (iounit, "(T12,1(F20.12,','),F20.12,' ]')") hessian(nc, j:nc)
ELSEIF (nc == j) THEN
WRITE (iounit, "(T12,F20.12,' ]')") hessian(nc, j:nc)
END IF
END DO
nc = 3*(i - 1) + 3
WRITE (iounit, "(T6,'Z: [')")
DO j = 1, nc, 5
IF (nc > j + 4) THEN
WRITE (iounit, "(T12,5(F20.12,','))") hessian(nc, j:j + 4)
ELSEIF (nc == j + 4) THEN
WRITE (iounit, "(T12,4(F20.12,','),F20.12,' ]')") hessian(nc, j:nc)
ELSEIF (nc == j + 3) THEN
WRITE (iounit, "(T12,3(F20.12,','),F20.12,' ]')") hessian(nc, j:nc)
ELSEIF (nc == j + 2) THEN
WRITE (iounit, "(T12,2(F20.12,','),F20.12,' ]')") hessian(nc, j:nc)
ELSEIF (nc == j + 1) THEN
WRITE (iounit, "(T12,1(F20.12,','),F20.12,' ]')") hessian(nc, j:nc)
ELSEIF (nc == j) THEN
WRITE (iounit, "(T12,F20.12,' ]')") hessian(nc, j:nc)
END IF
END DO
END DO
END IF
END IF

END SUBROUTINE write_scine

END MODULE scine_utils
@@ -0,0 +1,53 @@
&FORCE_EVAL
METHOD Quickstep
&PRINT
&SCINE
&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
A 6.0 0.0 0.0
B 0.0 6.0 1.0
C 1.0 0.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
@@ -4,4 +4,7 @@ C2H4_GRRM.inp 1 1.0E-12 -
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
#
C2H4_SCINE.inp 1 1.0E-12 -13.137211552379020
h2o_SCINE.inp 8 8.0E-06 1465.672
#EOF
@@ -0,0 +1,59 @@
&FORCE_EVAL
METHOD Quickstep
&PRINT
&SCINE
FILENAME =vib.scine
&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


0 comments on commit 688071a

Please sign in to comment.
You can’t perform that action at this time.