comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-29 06:35:28 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

New version:

Browse source 
Contains

1) fit modules
2) polarization modules

Includes many bug fixes to SPME.  There still seems to be a problem
with pw_integral_ab in regards to SPME/PME.  This version will
work for SPME not for PME.  Needs to be examined.


svn-origin-rev: 19
This commit is contained in:
Chris Mundy 2001-06-29 06:58:13 +00:00
parent 4ba9c53263
commit 638dbd3c61
64 changed files with 8350 additions and 1019 deletions
Download patch file Download diff file Expand all files Collapse all files

8
src/Makefile
View file

@ -1,9 +1,9 @@
.SUFFIXES: .F .d .o
CP2KHOME = $(HOME)/CP2K/
CP2KHOME = $(HOME)/CP2K
ARCH = $(shell $(CP2KHOME)/tools/get_arch_code)
VERSION = sopt
SMAKE = make -r
SMAKE = gmake -r
PMAKE = $(SMAKE) -j 1
LIB1 = cp2klib
@ -81,7 +81,7 @@ $(PROG): $(OBJECTS)
$(CPP) $(CPPFLAGS) $< $*.f
$(FORPAR) $*.f
$(PERL) $(SFMAKEDEPEND) $*.d $*.f
@rm -f $*.f $*.d.old
@rm -f $*.d.old
parallel_include.o: parallel_include.F
$(FC_fixed) -c $(FCFLAGS) $<
@ -89,7 +89,7 @@ parallel_include.o: parallel_include.F
parallel_include.d: parallel_include.F
$(CPP) $(CPPFLAGS) $< $*.f
$(PERL) $(SFMAKEDEPEND) $*.d $*.f
@rm -f $*.f $*.d.old
@rm -f $*.d.old
clean:
rm -f *.o *.mod F*.f

23
src/OBJECTDEFS
View file

@ -6,6 +6,7 @@ OBJECTS_GENERIC =\
ai_overlap_ppl.o\
ai_verfc.o\
amoeba.o\
ao_types.o\
atom_module.o\
atomic_kinds.o\
atoms.o\
@ -16,6 +17,7 @@ OBJECTS_GENERIC =\
brillouin.o\
cell_parameters.o\
cntl_input.o\
coefficient_lists.o\
coefficient_types.o\
coefficients.o\
collocate_density.o\
@ -42,6 +44,7 @@ OBJECTS_GENERIC =\
ewald_parameters_types.o\
ewalds.o\
external_potential_types.o\
extended_system_types.o\
fermi.o\
fft_tools.o\
fftessl_lib.o\
@ -53,9 +56,16 @@ OBJECTS_GENERIC =\
fist_force_numer.o\
fist_intra_force.o\
fist_nonbond_force.o\
fit_force_control.o\
fit_input.o\
fit_module.o\
fit_parameters_types.o\
fit_types.o\
fit_utilities.o\
force_control.o\
force_fields.o\
functionals.o\
function_amoeba.o\
gamma.o\
global_types.o\
greens_fn.o\
@ -64,7 +74,7 @@ OBJECTS_GENERIC =\
initialize_extended_types.o\
initialize_molecule_types.o\
initialize_particle_types.o\
initialize_pimd_types.o \
initialize_pimd_types.o\
input_types.o\
input_utilities.o\
integrate_potential.o\
@ -77,9 +87,11 @@ OBJECTS_GENERIC =\
kpoints.o\
lapack.o\
library_tests.o\
linear_systems.o\
linklist_cell_list.o\
linklist_cell_types.o\
linklist_control.o\
linklist_overlap_list.o\
linklist_types.o\
linklist_utilities.o\
linklist_verlet_list.o\
@ -117,6 +129,11 @@ OBJECTS_GENERIC =\
periodic_table.o\
physcon.o\
pme.o\
pol_coefs.o\
pol_electrostatics.o\
pol_force.o\
pol_grids.o\
pol_setup.o\
potential_types.o\
potentials.o\
print_keys.o\
@ -132,8 +149,8 @@ OBJECTS_GENERIC =\
radial_poisson.o\
radial_solver.o\
radial_util.o\
read_pimd.o\
readin.o\
read_pimd.o\
realspace_grid_types.o\
simulation_cell.o\
slater_koster.o\
@ -156,7 +173,7 @@ OBJECTS_GENERIC =\
termination.o\
timesl.o\
timings.o\
transformations.o \
transformations.o\
unit.o\
util.o\
wf_types.o\

77
src/amoeba.F
View file

@ -17,10 +17,20 @@
!! SOURCE
!******************************************************************************
MODULE amoeba
MODULE amoeba
USE kinds, ONLY : dbl
USE nrutil, ONLY : assert_eq, imaxloc, iminloc, nrerror, swap
USE structure_types, ONLY : structure_type
USE pair_potential, ONLY : potentialparm_type
USE particle_types, ONLY : particle_prop_type
USE md, ONLY : simulation_parameters_type, thermodynamic_type
USE fit_parameters_types, ONLY : fit_parameter_type, linklist_parameters_type, &
fit_info_type
USE ewald_parameters_types, ONLY : ewald_parameters_type
USE global_types, ONLY : global_environment_type
USE simulation_cell, ONLY : cell_type
IMPLICIT NONE
@ -49,29 +59,68 @@ CONTAINS
! NUMERICAL RECIPIES SUBROUTINE AMOEBA ADAPTED FOR USE IN CP2K !
!-----------------------------------------------------------------------------!
SUBROUTINE amoeba_evaluate ( p, y, ftol, rtol, func, itmax )
SUBROUTINE amoeba_evaluate ( p, y, rtol, struc, thermo, potparm, simpar, &
fit_param, box_ref, finfo, param, pstat, ewald_param, globenv, func )
IMPLICIT NONE
! Arguments
REAL (dbl), DIMENSION (:,:), INTENT (INOUT) :: p
REAL (dbl), DIMENSION (:), INTENT (INOUT) :: y
REAL (dbl), INTENT (IN) :: ftol
REAL (dbl), INTENT ( OUT ) :: rtol
INTEGER, INTENT ( IN ) :: itmax
TYPE (structure_type), INTENT (INOUT) :: struc
TYPE (potentialparm_type), DIMENSION (:,:), INTENT (INOUT) :: potparm
TYPE (thermodynamic_type), INTENT (INOUT) :: thermo
TYPE (simulation_parameters_type), INTENT (IN) :: simpar
TYPE (fit_parameter_type), INTENT(INOUT) :: fit_param
TYPE (cell_type) , POINTER, DIMENSION (:) :: box_ref
TYPE (fit_info_type), DIMENSION (:,:), POINTER :: finfo
TYPE (linklist_parameters_type ), INTENT(INOUT), DIMENSION ( :, : ) :: param
TYPE (particle_prop_type ), DIMENSION ( : ), POINTER :: pstat
TYPE ( ewald_parameters_type ), INTENT ( INOUT ) :: ewald_param
TYPE ( global_environment_type), INTENT ( IN ) :: globenv
! Interface
INTERFACE
FUNCTION func ( x )
USE kinds, ONLY : dbl
IMPLICIT NONE
REAL ( dbl ), DIMENSION ( : ), INTENT ( IN ) :: x
REAL ( dbl ) :: func
FUNCTION func ( x ,struc,potparm,thermo,simpar,fit_param,box_ref,&
finfo,param,pstat,ewald_param, globenv)
USE kinds, ONLY : dbl
USE structure_types, ONLY : structure_type
USE pair_potential, ONLY : potentialparm_type
USE particle_types, ONLY : particle_prop_type
USE md, ONLY : simulation_parameters_type, thermodynamic_type
USE fit_parameters_types, ONLY : fit_parameter_type, linklist_parameters_type, &
fit_info_type
USE ewald_parameters_types, ONLY : ewald_parameters_type
USE global_types, ONLY : global_environment_type
USE simulation_cell, ONLY : cell_type
IMPLICIT NONE
TYPE (structure_type), INTENT (INOUT) :: struc
TYPE (potentialparm_type), DIMENSION (:,:), INTENT (INOUT) :: potparm
TYPE (thermodynamic_type), INTENT (INOUT) :: thermo
TYPE (simulation_parameters_type), INTENT (IN) :: simpar
REAL (dbl), INTENT (IN), DIMENSION (:) :: x
TYPE (fit_parameter_type), INTENT(INOUT) :: fit_param
TYPE (cell_type) , POINTER, DIMENSION (:) :: box_ref
TYPE (fit_info_type), DIMENSION (:,:), POINTER :: finfo
TYPE (linklist_parameters_type ), INTENT(INOUT), DIMENSION ( :, : ) :: param
TYPE (particle_prop_type ), DIMENSION ( : ), POINTER :: pstat
TYPE ( ewald_parameters_type ), INTENT ( INOUT ) :: ewald_param
TYPE ( global_environment_type), INTENT ( IN ) :: globenv
REAL ( dbl ) :: func
END FUNCTION func
END INTERFACE
INTEGER :: ihi, ndim
REAL (dbl), DIMENSION (size(p,2)) :: psum
INTEGER :: ihi, ndim, itmax
REAL (dbl), DIMENSION (size(p,2)) :: psum
REAL (dbl) :: ftol
ftol = fit_param % max_tol
itmax = fit_param % iwrite
CALL amoeba_private()
@ -114,7 +163,8 @@ CONTAINS
IF (ytry>=ysave) THEN
p(:,:) = 0.5_dbl*(p(:,:)+spread(p(ilo,:),1,size(p,1)))
DO i = 1, ndim + 1
IF (i/=ilo) y(i) = func(p(i,:))
IF (i/=ilo) y(i) = func(p(i,:),struc,potparm,thermo,simpar,fit_param,&
box_ref,finfo,param, pstat, ewald_param, globenv)
END DO
iter = iter + ndim
psum(:) = sum(p(:,:),dim=1)
@ -134,7 +184,8 @@ CONTAINS
fac1 = (1.0_dbl-fac)/ndim
fac2 = fac1 - fac
ptry(:) = psum(:)*fac1 - p(ihi,:)*fac2
ytry = func(ptry)
ytry = func(ptry,struc,potparm,thermo,simpar,fit_param,&
box_ref,finfo,param, pstat, ewald_param, globenv)
IF (ytry<y(ihi)) THEN
y(ihi) = ytry
psum(:) = psum(:) - p(ihi,:) + ptry(:)

229
src/ao_types.F Normal file
View file

@ -0,0 +1,229 @@
!-----------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!***** cp2k/ao_types [1.0] *
!!
!! NAME
!! ao_types
!!
!! FUNCTION
!!
!! AUTHOR
!! CJM
!!
!! MODIFICATION HISTORY
!!
!! NOTES
!!
!!
!! SOURCE
!******************************************************************************
MODULE ao_types
USE atomic_kinds, ONLY : kind_info_type
USE kinds, ONLY : dbl
USE termination, ONLY : stop_memory, stop_program
USE timings, ONLY : timeset, timestop
IMPLICIT NONE
PRIVATE
PUBLIC :: ao_type
PUBLIC :: ao_allocate, ao_deallocate
PUBLIC :: ao_zero
TYPE ao_type
REAL ( dbl ), DIMENSION ( : ), POINTER :: cr
REAL ( dbl ), DIMENSION ( : ), POINTER :: mass
INTEGER, DIMENSION ( : ), POINTER :: coef_to_basis_set
INTEGER, DIMENSION ( : ), POINTER :: coef_to_part
INTEGER, DIMENSION ( : ), POINTER :: grid_id
TYPE ( kind_info_type ), DIMENSION ( : ), POINTER :: kind_info
END TYPE ao_type
!!*****
!******************************************************************************
CONTAINS
!******************************************************************************
!!****** ao_types/ao_allocate [1.0] *
!!
!! NAME
!! ao_allocate
!!
!! FUNCTION
!! Allocate a atomic oribital type
!!
!! AUTHOR
!! CJM
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
SUBROUTINE ao_allocate ( ao, ncoefs, mass_flag, grid_flag )
IMPLICIT NONE
! Arguments
TYPE ( ao_type ), INTENT ( INOUT ), TARGET :: ao
INTEGER, INTENT ( IN ) :: ncoefs
LOGICAL, INTENT ( IN ), OPTIONAL :: mass_flag, grid_flag
! Locals
INTEGER :: allocstat
LOGICAL :: local_mass_flag, local_grid_flag
!------------------------------------------------------------------------------
IF ( PRESENT ( mass_flag ) ) THEN
local_mass_flag=mass_flag
ELSE
local_mass_flag=.FALSE.
ENDIF
IF ( PRESENT ( grid_flag ) ) THEN
local_grid_flag=grid_flag
ELSE
local_grid_flag=.FALSE.
ENDIF
NULLIFY ( ao % cr, ao % coef_to_basis_set, ao % coef_to_part )
NULLIFY ( ao % mass )
NULLIFY ( ao % grid_id )
ALLOCATE ( ao % cr ( ncoefs ), STAT = allocstat )
IF ( allocstat /= 0 ) CALL stop_memory ( "ao_allocate", &
"ao % cr", ncoefs )
ALLOCATE ( ao % coef_to_basis_set ( ncoefs ), STAT = allocstat )
IF ( allocstat /= 0 ) CALL stop_memory ( "ao_allocate", &
"ao % coef_to_basis_set", ncoefs )
ALLOCATE ( ao % coef_to_part ( ncoefs ), STAT = allocstat )
IF ( allocstat /= 0 ) CALL stop_memory ( "ao_allocate", &
"ao % coef_to_part", ncoefs )
IF ( local_mass_flag ) THEN
ALLOCATE ( ao % mass ( ncoefs ), STAT = allocstat )
IF ( allocstat /= 0 ) CALL stop_memory ( "ao_allocate", &
"ao % mass", ncoefs )
ENDIF
IF ( local_grid_flag ) THEN
ALLOCATE ( ao % grid_id ( ncoefs ), STAT = allocstat )
IF ( allocstat /= 0 ) CALL stop_memory ( "ao_allocate", &
"ao % grid_id", ncoefs )
ENDIF
END SUBROUTINE ao_allocate
!!*****
!******************************************************************************
!!****** ao_types/ao_deallocate [1.0] *
!!
!! NAME
!! ao_deallocate
!!
!! FUNCTION
!! Deallocate a atomic orbital type variable
!!
!! AUTHOR
!! CJM
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
SUBROUTINE ao_deallocate ( ao )
IMPLICIT NONE
! Arguments
TYPE ( ao_type ), INTENT ( INOUT ) :: ao
! Locals
INTEGER :: allocstat
!------------------------------------------------------------------------------
DEALLOCATE ( ao % cr, STAT = allocstat )
IF ( allocstat /= 0 ) CALL stop_memory ( "ao_deallocate", &
"ao % cr", SIZE ( ao % cr ) )
DEALLOCATE ( ao % coef_to_basis_set, STAT = allocstat )
IF ( allocstat /= 0 ) CALL stop_memory ( "ao_deallocate", &
"ao % coef_to_basis_set", SIZE ( ao % coef_to_basis_set ) )
DEALLOCATE ( ao % coef_to_part, STAT = allocstat )
IF ( allocstat /= 0 ) CALL stop_memory ( "ao_deallocate", &
"ao % coef_to_part", SIZE ( ao % coef_to_part ) )
IF ( ASSOCIATED ( ao % mass )) THEN
DEALLOCATE ( ao % mass, STAT = allocstat )
IF ( allocstat /= 0 ) CALL stop_memory ( "ao_deallocate", &
"ao % mass", SIZE ( ao % mass ) )
ENDIF
IF ( ASSOCIATED ( ao % grid_id )) THEN
DEALLOCATE ( ao % mass, STAT = allocstat )
IF ( allocstat /= 0 ) CALL stop_memory ( "ao_deallocate", &
"ao % grid_id", SIZE ( ao % grid_id ) )
ENDIF
END SUBROUTINE ao_deallocate
!!*****
!******************************************************************************
!!****** ao_types/ao_zero [1.0] *
!!
!! NAME
!! ao_zero
!!
!! FUNCTION
!! Set values of a ao type to zero
!!
!! AUTHOR
!! CJM
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
SUBROUTINE ao_zero ( ao )
IMPLICIT NONE
! Arguments
TYPE ( ao_type ), INTENT ( INOUT ) :: ao
! Local
INTEGER :: ns, handle
REAL ( dbl ) :: zr
!------------------------------------------------------------------------------
CALL timeset ( 'AO_ZERO', 'I', 'Mzero', handle )
ns = SIZE ( ao % cr)
ao % cr = 0.0_dbl
ao % coef_to_basis_set = 0
IF ( ASSOCIATED ( ao % mass ) ) ao % mass = 0.0_dbl
IF ( ASSOCIATED ( ao % grid_id ) ) ao % grid_id = 0
zr = REAL ( ns, dbl ) * 1.e-6_dbl
CALL timestop ( zr, handle )
END SUBROUTINE ao_zero
!!*****
!******************************************************************************
END MODULE ao_types
!******************************************************************************

1
src/atomic_kinds.F
View file

@ -76,6 +76,7 @@ MODULE atomic_kinds
core_charge_radius,mass,&
prj_paw_radius
INTEGER :: natom,z,zeff
INTEGER :: number_of_grid_points
INTEGER, DIMENSION(:), POINTER :: atom_list,elec_conf
END TYPE kind_info_type

33
src/atoms_input.F
View file

@ -14,6 +14,7 @@
!! MODIFICATION HISTORY
!! cjm, Feb-20-2001 : added all the extended variables to
!! system_type
!! gt, Apr-05-2001 : atype pointers nullified
!!
!! SOURCE
!******************************************************************************
@ -144,13 +145,13 @@ SUBROUTINE read_coord_vel ( atype, filen, globenv )
CALL mp_bcast ( atype % ptype, source, allgrp )
CALL mp_bcast ( atype % rtype, source, allgrp )
CALL mp_bcast ( atype % n, source, allgrp )
IF ( .NOT. ASSOCIATED ( atype % c ) ) THEN
IF ( .NOT. globenv % ionode ) THEN
ALLOCATE ( atype % c(1:3,1:atype % n ), STAT = ios )
IF ( ios /= 0 ) CALL stop_memory ( 'ATOM', 'atype%c', 3 * atype % n )
END IF
CALL mp_bcast ( atype % c, source, allgrp )
IF ( atype % rtype == 'POSVEL' ) THEN
IF ( .NOT. ASSOCIATED ( atype % v ) ) THEN
IF ( .NOT. globenv % ionode ) THEN
ALLOCATE (atype % v(1:3,1:atype % n),STAT=ios)
IF ( ios /= 0 ) CALL stop_memory ( 'ATOM', 'atype%v', 3 * atype % n )
END IF
@ -257,13 +258,9 @@ SUBROUTINE read_file ( filen, atype )
OPEN ( iunit, file = filen )
READ ( iunit, * ) atype % n
IF (.NOT. ASSOCIATED (atype % c) ) THEN
ALLOCATE ( atype % c ( 1:3, 1:atype % n ), STAT = ios )
IF ( ios /= 0 ) CALL stop_memory ( 'ATOM', 'atype%c', 3 * atype % n )
NULLIFY ( atype % v )
ELSE IF (atype % n /= SIZE( atype % c, 2 ) ) THEN
CALL stop_program ( 'atoms_input','Inconsistent number of atoms' )
END IF
ALLOCATE ( atype % c ( 1:3, 1:atype % n ), STAT = ios )
IF ( ios /= 0 ) CALL stop_memory ( 'ATOM', 'atype%c', 3 * atype % n )
NULLIFY ( atype % v )
IF ( atype % rtype == 'INIT' ) THEN
DO i = 1, atype % n
@ -285,11 +282,9 @@ SUBROUTINE read_file ( filen, atype )
READ ( iunit, * ) atype % box_ref ( 2, 1:3 )
READ ( iunit, * ) atype % box_ref ( 3, 1:3 )
ELSE IF ( atype % rtype == 'POSVEL' ) THEN
IF (.NOT. ASSOCIATED (atype % v) ) THEN
ALLOCATE ( atype % v ( 1:3, 1:atype % n ), STAT = ios )
IF ( ios /= 0 ) &
CALL stop_memory ( 'ATOM', 'atype%v', 3 * atype % n )
END IF
ALLOCATE ( atype % v ( 1:3, 1:atype % n ), STAT = ios )
IF ( ios /= 0 ) &
CALL stop_memory ( 'ATOM', 'atype%v', 3 * atype % n )
DO i = 1, atype % n
READ ( iunit, * ) atype % c ( 1:3, i )
@ -305,11 +300,9 @@ SUBROUTINE read_file ( filen, atype )
READ ( iunit, * ) atype % v ( 1:3, i )
END DO
ELSE IF ( atype % rtype == 'ALL' ) THEN
IF (.NOT. ASSOCIATED (atype % v) ) THEN
ALLOCATE ( atype % v ( 1:3, 1:atype % n ), STAT = ios )
IF ( ios /= 0 ) &
CALL stop_memory ( 'ATOM', 'atype%v', 3 * atype % n )
END IF
ALLOCATE ( atype % v ( 1:3, 1:atype % n ), STAT = ios )
IF ( ios /= 0 ) &
CALL stop_memory ( 'ATOM', 'atype%v', 3 * atype % n )
DO i = 1, atype % n
READ ( iunit, * ) atype % c ( 1:3, i )
@ -404,8 +397,6 @@ SUBROUTINE read_file ( filen, atype )
CALL stop_program ( 'ATOM', 'this rtype not programmed' )
END IF
CLOSE ( iunit )
ELSE
CALL stop_program ( 'ATOM', 'No information on atoms found ' )

2
src/cntl_input.F
View file

@ -52,7 +52,7 @@ CONTAINS
!!I---------------------------------------------------------------------------I
!!I SECTION: &cntl ... &end I
!!I I
!!I simulation [md,pimd,debug] I
!!I simulation [md,pimd,debug,min,amoeba] I
!!I printlevel globenv%print_level I
!!I units [kelvin,atomic] I
!!I periodic [0,1][0,1][0,1] I

112
src/coefficient_lists.F Normal file
View file

@ -0,0 +1,112 @@
!-----------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/coefficient_list [1.0] *
!!
!! NAME
!! coefficient_lists
!!
!! FUNCTION
!!
!! AUTHOR
!! gloria
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE coefficient_lists
USE kinds, ONLY : dbl
USE particle_types, ONLY : particle_type
use ao_types, ONLY : ao_type
USE termination, ONLY : stop_memory
IMPLICIT NONE
PRIVATE
PUBLIC :: coef_list_type, coef_list_grid
TYPE coef_list_type
INTEGER :: coef_index
INTEGER :: particle_index
INTEGER :: grid_type_id
INTEGER :: grid_id
REAL ( dbl ), DIMENSION ( : ), POINTER :: r
REAL ( dbl ), POINTER :: c
TYPE ( coef_list_type ), POINTER :: next
END TYPE coef_list_type
!!*****
!******************************************************************************
CONTAINS
!******************************************************************************
SUBROUTINE coef_list_grid ( ao, clist, part, ncoef, grid_type_tag, &
grid_tag)
IMPLICIT NONE
! Arguments
TYPE ( ao_type ), INTENT ( IN ) :: ao
TYPE ( coef_list_type ), POINTER :: clist
TYPE ( particle_type ), intent (IN), DIMENSION (:), TARGET :: PART
INTEGER , INTENT ( OUT ) :: ncoef
INTEGER , INTENT ( IN ) :: grid_type_tag
INTEGER , INTENT ( IN ) :: grid_tag
! Locals
INTEGER :: ic, nc, isos, ipart
TYPE ( coef_list_type ), POINTER :: clist_loc
!------------------------------------------------------------------------------
nc = SIZE ( ao%cr )
! 'plist' will always point to the beginning of the list, 'plist_loc'
! to the current position
NULLIFY ( clist )
NULLIFY ( clist_loc )
ncoef = 0
DO ic = 1, nc
IF ( ao%grid_id ( ic ) == grid_type_tag ) THEN
IF ( .NOT. ASSOCIATED ( clist_loc ) ) THEN
ALLOCATE ( clist, STAT = isos )
IF ( isos /= 0 ) &
CALL stop_memory ( "get_atom_list", "plist", 0 )
clist_loc => clist
NULLIFY ( clist_loc % next )
ELSE
ALLOCATE ( clist_loc % next, STAT = isos )
IF ( isos /= 0 ) &
CALL stop_memory ( "get_atom_list", "next", 0 )
clist_loc => clist_loc % next
NULLIFY ( clist_loc % next )
END IF
ncoef = ncoef + 1
clist_loc % coef_index = ic
ipart = ao % coef_to_part ( ic )
clist_loc % particle_index = ipart
clist_loc % grid_type_id = grid_type_tag
clist_loc % grid_id = grid_tag
clist_loc % r => part ( ipart ) % r
clist_loc % c => ao % cr (ic)
END IF
END DO
END SUBROUTINE coef_list_grid
!******************************************************************************
END MODULE coefficient_lists

8
src/coefficient_types.F
View file

@ -21,6 +21,7 @@
MODULE coefficient_types
USE kinds, ONLY : dbl
USE ao_types, ONLY : ao_type
USE pw_grid_types, ONLY : pw_grid_type
USE pw_types, ONLY : pw_type, pw_zero, pw_copy, &
pw_add, pw_sumup, pw_subtract, pw_reduce, &
@ -43,6 +44,7 @@ MODULE coefficient_types
TYPE coeff_type
INTEGER :: in_use
TYPE ( pw_type ) :: pw
TYPE ( ao_type ) :: ao
END TYPE coeff_type
! Flags for the structure member 'in_use'
@ -344,7 +346,11 @@ SUBROUTINE coeff_transform_c1c2 ( c1, c2, debug )
!..check if we are dealing with plane waves and call the correct routine
IF ( c1 % in_use == PLANEWAVES .AND. c2 % in_use == PLANEWAVES ) THEN
CALL pw_fft_wrap ( c1 % pw, c2 % pw, debug )
IF ( PRESENT ( debug ) ) THEN
CALL pw_fft_wrap ( c1 % pw, c2 % pw, debug )
ELSE
CALL pw_fft_wrap ( c1 % pw, c2 % pw )
END IF
ELSE
CALL stop_program ( "coeff_transform_c1c2", "not implemented")
END IF

966
src/control_module.F
View file

File diff suppressed because it is too large Load diff

2
src/cp2k.F
View file

@ -78,7 +78,7 @@ PROGRAM cp2k
CALL stop_program ( "cp2k", "No valid program name given "&
//globenv % program_name )
CASE ( "FIST", "TBMD", "WAVE" )
CASE ( "FIST", "TBMD", "WAVE", "FIT", "POL" )
CALL control ( globenv )

13
src/dft_structure_control.F
View file

@ -29,7 +29,6 @@ MODULE dft_structure_control
USE dft_input, ONLY : read_dft_section
USE ewald_parameters_types, ONLY : ewald_parameters_type
USE global_types, ONLY : global_environment_type
USE greens_fn, ONLY : pw_green_fn_init
USE input_types, ONLY : setup_parameters_type
USE kinds, ONLY : dbl, YES, NO
USE linklist_control, ONLY : set_ll_parm
@ -47,7 +46,6 @@ MODULE dft_structure_control
USE potentials, ONLY : potential_update
USE pw_grids, ONLY : pw_grid_setup, pw_grid_construct
USE pw_grid_types, ONLY : pw_grid_type, HALFSPACE, FULLSPACE
USE pws, ONLY : init_pw_poisson_solver
USE pw_types, ONLY : PW_REALDATA3D => REALDATA3D
USE simulation_cell, ONLY : cell_type, get_hinv
USE termination, ONLY : stop_memory
@ -99,8 +97,6 @@ SUBROUTINE dft_force_control ( molecule, pnode, part, box, thermo, &
ALLOCATE ( tau ( 3, 1 ), fion ( 3, 1 ) ) !TMPTMPTMP
CALL init_pw_poisson_solver ( box % green )
DO
CALL dft_scf_loop ( psi, n, vloc, tau, fion, dft_control, kp, &
dftpar, wavepar )
@ -233,11 +229,6 @@ SUBROUTINE dft_initialise ( cell, psi, n, vloc, dft_control, kp, wavepar )
CALL wf_randomise ( psi_n, wavepar )
END DO
! the Greens function
! choose the set that is used to solve the Poisson equation in G space
densityset = 1
CALL pw_green_fn_init ( cell % green, pw_grid ( densityset ) )
END SUBROUTINE dft_initialise
@ -248,8 +239,6 @@ SUBROUTINE dft_finalise ( dft_control )
TYPE ( dft_control_type ) :: dft_control
!!!!! if the pw_grid is deallocated -> do not forget to get rid of the Greens function
END SUBROUTINE dft_finalise
!******************************************************************************
@ -353,5 +342,3 @@ END SUBROUTINE dft_scf_loop
!******************************************************************************
END MODULE dft_structure_control
!******************************************************************************

255
src/dg_rho0s.F
View file

@ -20,7 +20,7 @@
MODULE dg_rho0s
USE coefficient_types, ONLY : coeff_type, coeff_allocate, coeff_zero, &
PW_REALDATA3D
PW_REALDATA3D, PW_COMPLEXDATA3D
USE kinds, ONLY : dbl
USE mathconstants, ONLY : pi
USE pw_grid_types, ONLY : pw_grid_type
@ -38,6 +38,8 @@ MODULE dg_rho0s
TYPE dg_rho0_type
CHARACTER ( LEN = 12 ) :: type
INTEGER :: grid
INTEGER :: kind
REAL ( dbl ) :: cutoff_radius
REAL ( dbl ), DIMENSION ( : ), POINTER :: gcc
REAL ( dbl ), DIMENSION ( : ), POINTER :: zet
TYPE ( coeff_type ) :: density
@ -70,7 +72,16 @@ SUBROUTINE dg_rho0_setup ( dg_rho0, pw_grid )
CALL coeff_allocate ( dg_rho0 % density, pw_grid, PW_REALDATA3D )
CALL dg_rho0_pme_gauss ( dg_rho0 % density, dg_rho0 % zet ( 1 ) )
CASE ( "SPME" )
CALL stop_program ( 'dg_rho0_setup', 'type not implemented' )
CALL coeff_allocate ( dg_rho0 % density, pw_grid, PW_REALDATA3D )
CALL dg_rho0_pme_gauss ( dg_rho0 % density, dg_rho0 % zet ( 1 ) )
CASE ( "S" )
CALL coeff_allocate ( dg_rho0 % density, pw_grid, PW_REALDATA3D )
CALL dg_rho0_s_cgf ( dg_rho0 % density, dg_rho0 % zet , &
dg_rho0 % gcc)
CASE ( "PX", "PY", "PZ" )
CALL coeff_allocate ( dg_rho0 % density, pw_grid, PW_COMPLEXDATA3D )
CALL dg_rho0_p_cgf ( dg_rho0 % density, dg_rho0 % zet, &
dg_rho0 % gcc, dg_rho0 % type)
CASE DEFAULT
CALL stop_program ( 'dg_rho0_setup', 'type not implemented' )
END SELECT
@ -157,6 +168,246 @@ SUBROUTINE dg_rho0_pme_gauss ( dg_rho0, alpha )
END DO
END SUBROUTINE dg_rho0_pme_gauss
!******************************************************************************
SUBROUTINE dg_rho0_s_cgf ( dg_rho0, zet, gcc )
!******************************************************************************
! NAME
!! dg_rho0_s_cgf
!!
!! FUNCTION
!! puts a contracted gaussian function of s-type
!! in the g-space
!!
!! AUTHOR
!! gloria
!!
!! SOURCE
!******************************************************************************
IMPLICIT NONE
! Arguments
TYPE ( coeff_type ), INTENT ( INOUT ), TARGET :: dg_rho0
REAL ( dbl ), INTENT ( IN ), DIMENSION ( : ) :: zet
REAL ( dbl ), INTENT ( IN ), DIMENSION ( : ) :: gcc
! Locals
TYPE ( pw_grid_type ), POINTER :: pw_grid
INTEGER :: l0, m0, n0
INTEGER, PARAMETER :: IMPOSSIBLE = 10000
REAL ( dbl ) :: const, e_gsq, charge
INTEGER :: lp, mp, np, ln, mn, nn, gpt, ncg, icg, coef
INTEGER, DIMENSION ( :, : ), POINTER :: bds
INTEGER, DIMENSION ( : ), POINTER :: ghat
REAL ( dbl ), DIMENSION ( :, :, : ), POINTER :: rho0
!------------------------------------------------------------------------------
ncg = SIZE ( zet ( : ) )
pw_grid => dg_rho0 % pw % pw_grid
bds => pw_grid % bounds
IF ( -bds ( 1, 1 ) == bds ( 2, 1 ) ) THEN
l0 = IMPOSSIBLE
ELSE
l0 = bds ( 1, 1 )
END IF
IF ( -bds ( 1, 2 ) == bds ( 2, 2 ) ) THEN
m0 = IMPOSSIBLE
ELSE
m0 = bds ( 1, 2 )
END IF
IF ( -bds ( 1, 3 ) == bds ( 2, 3 ) ) THEN
n0 = IMPOSSIBLE
ELSE
n0 = bds ( 1, 3 )
END IF
CALL coeff_zero ( dg_rho0 )
rho0 => dg_rho0 % pw % cr3d
charge = 0._dbl
DO icg = 1, ncg
const = 1.0_dbl / ( 4.0_dbl * zet ( icg ) )
coef = ( gcc ( icg ) * pi ** 1.5 ) / zet ( icg )**1.5
charge = charge + coef * 1._dbl
DO gpt = 1, pw_grid % ngpts_cut
ghat => pw_grid % g_hat ( :, gpt )
lp = pw_grid % mapl % pos ( ghat ( 1 ) )
ln = pw_grid % mapl % neg ( ghat ( 1 ) )
mp = pw_grid % mapm % pos ( ghat ( 2 ) )
mn = pw_grid % mapm % neg ( ghat ( 2 ) )
np = pw_grid % mapn % pos ( ghat ( 3 ) )
nn = pw_grid % mapn % neg ( ghat ( 3 ) )
e_gsq = EXP ( -const * pw_grid % gsq ( gpt ) ) / pw_grid % vol
!*apsi
lp = lp + bds ( 1, 1 )
mp = mp + bds ( 1, 2 )
np = np + bds ( 1, 3 )
ln = ln + bds ( 1, 1 )
mn = mn + bds ( 1, 2 )
nn = nn + bds ( 1, 3 )
rho0 ( lp, mp, np ) = rho0 ( lp, mp, np ) + e_gsq * coef
rho0 ( ln, mn, nn ) = rho0 ( ln, mn, nn ) + e_gsq * coef
IF ( ghat ( 1 ) == l0 .OR. ghat ( 2 ) == m0 .OR. ghat ( 3 ) == n0 ) THEN
rho0 ( lp, mp, np ) = 0.0_dbl
rho0 ( ln, mn, nn ) = 0.0_dbl
END IF
END DO
END DO
END SUBROUTINE dg_rho0_s_cgf
!******************************************************************************
SUBROUTINE dg_rho0_p_cgf ( dg_rho0, zet, gcc, type )
!******************************************************************************
! NAME
!! dg_rho0_p_cgf
!!
!! FUNCTION
!! puts a contracted gaussian function of p-type
!! in the g-space. The p is imaginary, so
!! the p-function has then to be multiplied by (0.0, -1.0)
!! in order to get the true p.
!!
!! AUTHOR
!! gloria
!!
!! SOURCE
!******************************************************************************
IMPLICIT NONE
! Arguments
TYPE ( coeff_type ), INTENT ( INOUT ), TARGET :: dg_rho0
REAL ( dbl ), INTENT ( IN ), DIMENSION (:) :: zet
REAL ( dbl ), INTENT ( IN ), DIMENSION (:) :: gcc
CHARACTER ( LEN = 12 ), INTENT ( IN ) :: type
! Locals
TYPE ( pw_grid_type ), POINTER :: pw_grid
INTEGER :: l0, m0, n0
INTEGER, PARAMETER :: IMPOSSIBLE = 10000
REAL ( dbl ) :: charge, const, e_gsq, p_on_grid
INTEGER :: lp, mp, np, ln, mn, nn, gpt, ncg, icg
INTEGER, DIMENSION ( :, : ), POINTER :: bds
INTEGER, DIMENSION ( : ), POINTER :: ghat
COMPLEX ( dbl ), DIMENSION ( :, :, : ), POINTER :: rho0
COMPLEX ( dbl ) :: minus_i
REAL ( dbl ) :: coef
!------------------------------------------------------------------------------
minus_i = ( 0._dbl, -1.0_dbl )
ncg = size(zet(:))
pw_grid => dg_rho0 % pw % pw_grid
bds => pw_grid % bounds
IF ( -bds ( 1, 1 ) == bds ( 2, 1 ) ) THEN
l0 = IMPOSSIBLE
ELSE
l0 = bds ( 1, 1 )
END IF
IF ( -bds ( 1, 2 ) == bds ( 2, 2 ) ) THEN
m0 = IMPOSSIBLE
ELSE
m0 = bds ( 1, 2 )
END IF
IF ( -bds ( 1, 3 ) == bds ( 2, 3 ) ) THEN
n0 = IMPOSSIBLE
ELSE
n0 = bds ( 1, 3 )
END IF
CALL coeff_zero ( dg_rho0 )
rho0 => dg_rho0 % pw % cc3d
charge = 0._dbl
DO icg = 1, ncg
IF ( zet (icg) == 0 ) CYCLE
const = 1.0_dbl / ( 4.0_dbl * zet ( icg ) )
coef = ( gcc ( icg ) * pi**1.5 ) / ( 2.0_dbl * zet ( icg )**2.5 )
charge = charge + coef * 1._dbl
DO gpt = 1, pw_grid % ngpts_cut
ghat => pw_grid % g_hat ( :, gpt )
lp = pw_grid % mapl % pos ( ghat ( 1 ) )
ln = pw_grid % mapl % neg ( ghat ( 1 ) )
mp = pw_grid % mapm % pos ( ghat ( 2 ) )
mn = pw_grid % mapm % neg ( ghat ( 2 ) )
np = pw_grid % mapn % pos ( ghat ( 3 ) )
nn = pw_grid % mapn % neg ( ghat ( 3 ) )
e_gsq = EXP ( -const * pw_grid % gsq ( gpt ) ) / pw_grid % vol
!*apsi
lp = lp + bds ( 1, 1 )
mp = mp + bds ( 1, 2 )
np = np + bds ( 1, 3 )
ln = ln + bds ( 1, 1 )
mn = mn + bds ( 1, 2 )
nn = nn + bds ( 1, 3 )
SELECT CASE (type)
CASE ("PX")
p_on_grid = pw_grid % g (1, gpt) * e_gsq * coef
CASE ("PY")
p_on_grid = pw_grid % g (2, gpt) * e_gsq * coef
CASE ("PZ")
p_on_grid = pw_grid % g (3, gpt) * e_gsq * coef
END SELECT
rho0 ( lp, mp, np ) = rho0 ( lp, mp, np ) + p_on_grid * minus_i
rho0 ( ln, mn, nn ) = rho0 ( ln, mn, nn ) - p_on_grid * minus_i
IF ( ghat ( 1 ) == l0 .OR. ghat ( 2 ) == m0 .OR. ghat ( 3 ) == n0 ) THEN
rho0 ( lp, mp, np ) = 0.0_dbl
rho0 ( ln, mn, nn ) = 0.0_dbl
END IF
END DO
END DO
END SUBROUTINE dg_rho0_p_cgf
!****************************************************************************
!******************************************************************************

3
src/dg_types.F
View file

@ -19,6 +19,7 @@ MODULE dg_types
USE dg_rho0s, ONLY : dg_rho0_type
USE particle_lists, ONLY : particle_list_type
USE coefficient_lists, ONLY : coef_list_type
IMPLICIT NONE
@ -26,8 +27,10 @@ MODULE dg_types
TYPE dg_type
INTEGER :: grid_index
INTEGER :: nparts
INTEGER :: ncoefs
TYPE ( dg_rho0_type ) :: dg_rho0
TYPE ( particle_list_type ), POINTER :: plist
TYPE ( coef_list_type ), POINTER :: clist
END TYPE dg_type
PRIVATE

416
src/dgs.F
View file

@ -20,7 +20,8 @@
MODULE dgs
USE coefficient_types, ONLY : coeff_type
USE coefficient_types, ONLY : coeff_type, PW_REALDATA1D, &
PW_COMPLEXDATA1D, PW_REALDATA3D, PW_COMPLEXDATA3D
USE fft_tools, ONLY : fft_radix_operations, fft3d, FWFFT, BWFFT, &
FFT_RADIX_ALLOWED, FFT_RADIX_DISALLOWED, FFT_RADIX_NEXT, &
FFT_RADIX_CLOSEST
@ -40,15 +41,20 @@ MODULE dgs
PRIVATE
PUBLIC :: dg_get_patch, dg_get_delta
PUBLIC :: dg_grid_setup, dg_find_max_cutoff, dg_sum_patch, &
dg_sum_patch_force, dg_get_strucfac, dg_grid_change
PUBLIC :: dg_pol_grid_setup, dg_pme_grid_setup, dg_find_max_cutoff, &
dg_sum_patch, dg_sum_patch_force_3d, dg_sum_patch_force_1d, &
dg_get_strucfac, dg_grid_change, dg_find_cutoff
INTERFACE dg_sum_patch
MODULE PROCEDURE dg_sum_patch_coef, dg_sum_patch_arr
END INTERFACE
INTERFACE dg_sum_patch_force
MODULE PROCEDURE dg_sum_patch_force_coef, dg_sum_patch_force_arr
INTERFACE dg_sum_patch_force_3d
MODULE PROCEDURE dg_sum_patch_force_coef_3d, dg_sum_patch_force_arr_3d
END INTERFACE
INTERFACE dg_sum_patch_force_1d
MODULE PROCEDURE dg_sum_patch_force_coef_1d, dg_sum_patch_force_arr_1d
END INTERFACE
INTERFACE dg_get_patch
@ -59,8 +65,12 @@ MODULE dgs
MODULE PROCEDURE dg_add_patch_simple, dg_add_patch_folded
END INTERFACE
INTERFACE dg_int_patch
MODULE PROCEDURE dg_int_patch_simple, dg_int_patch_folded
INTERFACE dg_int_patch_3d
MODULE PROCEDURE dg_int_patch_simple_3d, dg_int_patch_folded_3d
END INTERFACE
INTERFACE dg_int_patch_1d
MODULE PROCEDURE dg_int_patch_simple_1d, dg_int_patch_folded_1d
END INTERFACE
!!*****
@ -70,7 +80,7 @@ CONTAINS
!******************************************************************************
SUBROUTINE dg_grid_setup ( box_b, npts_s, epsilon, alpha, grid_s, grid_b, &
SUBROUTINE dg_pme_grid_setup ( box_b, npts_s, epsilon, alpha, grid_s, grid_b, &
dg_gaussian_type, dg_glopar )
IMPLICIT NONE
@ -129,9 +139,112 @@ SUBROUTINE dg_grid_setup ( box_b, npts_s, epsilon, alpha, grid_s, grid_b, &
CALL pw_grid_setup ( box_s, grid_s, symm_usage=.FALSE. )
END SUBROUTINE dg_grid_setup
END SUBROUTINE dg_pme_grid_setup
!******************************************************************************
SUBROUTINE dg_pol_grid_setup ( box_b, npts_s, cutoff_radius, &
grid_s, grid_b, dg_glopar, grid_ref )
IMPLICIT NONE
! Arguments
INTEGER, DIMENSION ( : ), INTENT ( IN ) :: npts_s
TYPE ( cell_type ), INTENT ( IN ) :: box_b
REAL ( dbl ) , INTENT ( IN ) :: cutoff_radius
TYPE ( global_environment_type ), INTENT ( IN ) :: dg_glopar
TYPE ( pw_grid_type ), INTENT ( INOUT ) :: grid_s, grid_b
TYPE ( pw_grid_type ), INTENT ( IN ), OPTIONAL :: grid_ref
! Locals
REAL ( dbl ) :: cutoff, ecut
TYPE ( cell_type ) :: unit_box, box_s
!------------------------------------------------------------------------------
CALL dg_find_cutoff ( box_b, npts_s, cutoff_radius, grid_s, grid_b, cutoff )
! If necessary set up grid with respect to the reference grid
ecut = 0.5_dbl * cutoff * cutoff
IF ( PRESENT ( grid_ref ) ) THEN
CALL pw_grid_setup ( box_b, grid_b, ecut, dg_glopar % group, &
ref_grid = grid_ref )
ELSE
CALL pw_grid_setup ( box_b, grid_b, ecut, dg_glopar % group )
ENDIF
CALL dg_find_basis ( grid_b % npts, box_b, unit_box )
CALL dg_set_cell ( grid_s % npts, unit_box, box_s )
CALL pw_grid_setup ( box_s, grid_s )
END SUBROUTINE dg_pol_grid_setup
!******************************************************************************
SUBROUTINE dg_find_cutoff ( box_b, npts_s, cutoff_radius, grid_s, &
grid_b, cutoff )
IMPLICIT NONE
! Arguments
INTEGER, DIMENSION ( : ), INTENT ( IN ) :: npts_s
TYPE ( cell_type ), INTENT ( IN ) :: box_b
TYPE ( pw_grid_type ), INTENT ( INOUT ) :: grid_s, grid_b
REAL ( dbl ) , INTENT ( IN ) :: cutoff_radius
REAL ( dbl ) , INTENT ( OUT ) :: cutoff
! Locals
INTEGER :: nout ( 3 )
REAL ( dbl ) :: dr ( 3 )
REAL ( dbl ) :: cell_lengths ( 3 )
CALL fft_radix_operations ( npts_s ( 1 ), nout ( 1 ), &
operation = FFT_RADIX_NEXT )
CALL fft_radix_operations ( npts_s ( 1 ), nout ( 2 ), &
operation = FFT_RADIX_NEXT )
CALL fft_radix_operations ( npts_s ( 1 ), nout ( 3 ), &
operation = FFT_RADIX_NEXT )
CALL get_cell_param ( box_b, cell_lengths )
CALL dg_pol_get_spacing ( nout, cutoff_radius, dr )
CALL dg_find_radix ( dr, cell_lengths, grid_b % npts )
grid_b % bounds ( 1, : ) = - grid_b % npts / 2
grid_b % bounds ( 2, : ) = + ( grid_b % npts - 1 ) / 2
grid_b % grid_span = HALFSPACE
grid_s % bounds ( 1, : ) = -nout ( : ) / 2
grid_s % bounds ( 2, : ) = ( +nout ( : ) - 1 ) / 2
grid_s % grid_span = HALFSPACE
grid_s % npts = nout
CALL pw_find_cutoff ( grid_b % npts, box_b, cutoff )
END SUBROUTINE dg_find_cutoff
!******************************************************************************
SUBROUTINE dg_pol_get_spacing ( npts, cutoff_radius, dr )
IMPLICIT NONE
! Arguments
REAL ( dbl ), INTENT ( IN ) :: cutoff_radius
INTEGER, DIMENSION ( : ), INTENT ( IN ) :: npts
REAL ( dbl ), DIMENSION ( : ), INTENT ( OUT ) :: dr
! Locals
REAL ( dbl ) :: alphasq, norm
!------------------------------------------------------------------------------
dr ( : ) = cutoff_radius / ( REAL ( npts ( : ), dbl ) / 2.0_dbl )
END SUBROUTINE dg_pol_get_spacing
!******************************************************************************
SUBROUTINE dg_grid_change ( box_b, grid_b, grid_s )
@ -467,7 +580,6 @@ SUBROUTINE dg_sum_patch_coef ( rs, rhos, center )
END IF
END SUBROUTINE dg_sum_patch_coef
!******************************************************************************
SUBROUTINE dg_sum_patch_arr ( rs, rhos, center )
@ -487,7 +599,6 @@ SUBROUTINE dg_sum_patch_arr ( rs, rhos, center )
LOGICAL :: folded
!------------------------------------------------------------------------------
nb => rs % npts
ns ( 1 ) = SIZE ( rhos, 1 )
ns ( 2 ) = SIZE ( rhos, 2 )
@ -549,7 +660,7 @@ END SUBROUTINE dg_sum_patch_arr
!******************************************************************************
SUBROUTINE dg_sum_patch_force_arr ( drpot, rhos, center, force )
SUBROUTINE dg_sum_patch_force_arr_3d ( drpot, rhos, center, force )
IMPLICIT NONE
@ -618,22 +729,104 @@ SUBROUTINE dg_sum_patch_force_arr ( drpot, rhos, center, force )
END DO
IF ( folded ) THEN
CALL dg_int_patch ( drpot ( 1 ) % r, drpot ( 2 ) % r, drpot ( 3 ) % r, &
CALL dg_int_patch_3d ( drpot ( 1 ) % r, drpot ( 2 ) % r, drpot ( 3 ) % r, &
rhos, force, ns, &
drpot ( 1 ) % px, drpot ( 1 ) % py, drpot ( 1 ) % pz )
ELSE
nc ( 1 ) = drpot ( 1 ) % px ( 1 ) - 1
nc ( 2 ) = drpot ( 1 ) % py ( 1 ) - 1
nc ( 3 ) = drpot ( 1 ) % pz ( 1 ) - 1
CALL dg_int_patch ( drpot ( 1 ) % r, drpot ( 2 ) % r, drpot ( 3 ) % r, &
CALL dg_int_patch_3d ( drpot ( 1 ) % r, drpot ( 2 ) % r, drpot ( 3 ) % r, &
rhos, force, ns, nc )
END IF
END SUBROUTINE dg_sum_patch_force_arr
END SUBROUTINE dg_sum_patch_force_arr_3d
!******************************************************************************
SUBROUTINE dg_sum_patch_force_coef ( drpot, rhos, center, force )
SUBROUTINE dg_sum_patch_force_arr_1d ( drpot, rhos, center, force )
IMPLICIT NONE
! Arguments
TYPE ( realspace_grid_type ), TARGET, &
INTENT ( INOUT ) :: drpot
REAL ( dbl ), DIMENSION ( :,:,: ), INTENT ( IN ) :: rhos
INTEGER, DIMENSION ( 3 ), INTENT ( IN ) :: center
REAL ( dbl ), INTENT ( OUT ) :: force
! Locals
INTEGER, DIMENSION ( : ), POINTER :: nb
INTEGER, DIMENSION ( 3 ) :: ns, lb, ub
INTEGER, DIMENSION ( 3 ) :: nc
INTEGER :: i, ii, ia
LOGICAL :: folded
!------------------------------------------------------------------------------
nb => drpot % npts
ns ( 1 ) = SIZE ( rhos, 1 )
ns ( 2 ) = SIZE ( rhos, 2 )
ns ( 3 ) = SIZE ( rhos, 3 )
lb = -(ns-1) / 2
ub = lb + ns - 1
folded = .FALSE.
DO i = lb ( 1 ), ub ( 1 )
ia = i - lb ( 1 ) + 1
ii = center ( 1 ) + i - drpot % lb ( 1 )
IF ( ii < 0 ) THEN
drpot % px ( ia ) = ii + nb ( 1 ) + 1
folded = .TRUE.
ELSEIF ( ii >= nb ( 1 ) ) THEN
drpot % px ( ia ) = ii - nb ( 1 ) + 1
folded = .TRUE.
ELSE
drpot % px ( ia ) = ii + 1
ENDIF
END DO
DO i = lb ( 2 ), ub ( 2 )
ia = i - lb ( 2 ) + 1
ii = center ( 2 ) + i - drpot % lb ( 2 )
IF ( ii < 0 ) THEN
drpot % py ( ia ) = ii + nb ( 2 ) + 1
folded = .TRUE.
ELSEIF ( ii >= nb ( 2 ) ) THEN
drpot % py ( ia ) = ii - nb ( 2 ) + 1
folded = .TRUE.
ELSE
drpot % py ( ia ) = ii + 1
ENDIF
END DO
DO i = lb ( 3 ), ub ( 3 )
ia = i - lb ( 3 ) + 1
ii = center ( 3 ) + i - drpot % lb ( 3 )
IF ( ii < 0 ) THEN
drpot % pz ( ia ) = ii + nb ( 3 ) + 1
folded = .TRUE.
ELSEIF ( ii >= nb ( 3 ) ) THEN
drpot % pz ( ia ) = ii - nb ( 3 ) + 1
folded = .TRUE.
ELSE
drpot % pz ( ia ) = ii + 1
ENDIF
END DO
IF ( folded ) THEN
CALL dg_int_patch_1d ( drpot % r, rhos, force, ns, &
drpot % px, drpot % py, drpot % pz )
ELSE
nc ( 1 ) = drpot % px ( 1 ) - 1
nc ( 2 ) = drpot % py ( 1 ) - 1
nc ( 3 ) = drpot % pz ( 1 ) - 1
CALL dg_int_patch_1d ( drpot % r, rhos, force, ns, nc )
END IF
END SUBROUTINE dg_sum_patch_force_arr_1d
!******************************************************************************
SUBROUTINE dg_sum_patch_force_coef_3d ( drpot, rhos, center, force )
IMPLICIT NONE
@ -698,18 +891,96 @@ SUBROUTINE dg_sum_patch_force_coef ( drpot, rhos, center, force )
END DO
IF ( folded ) THEN
CALL dg_int_patch ( drpot ( 1 ) % r, drpot ( 2 ) % r, drpot ( 3 ) % r, &
CALL dg_int_patch_3d ( drpot ( 1 ) % r, drpot ( 2 ) % r, drpot ( 3 ) % r, &
rhos % pw % cr3d, force, ns, &
drpot ( 1 ) % px, drpot ( 1 ) % py, drpot ( 1 ) % pz )
ELSE
nc ( 1 ) = drpot ( 1 ) % px ( 1 ) - 1
nc ( 2 ) = drpot ( 1 ) % py ( 1 ) - 1
nc ( 3 ) = drpot ( 1 ) % pz ( 1 ) - 1
CALL dg_int_patch ( drpot ( 1 ) % r, drpot ( 2 ) % r, drpot ( 3 ) % r, &
CALL dg_int_patch_3d ( drpot ( 1 ) % r, drpot ( 2 ) % r, drpot ( 3 ) % r, &
rhos % pw % cr3d, force, ns, nc )
END IF
END SUBROUTINE dg_sum_patch_force_coef
END SUBROUTINE dg_sum_patch_force_coef_3d
!******************************************************************************
SUBROUTINE dg_sum_patch_force_coef_1d ( drpot, rhos, center, force )
IMPLICIT NONE
! Arguments
TYPE ( realspace_grid_type ), TARGET, INTENT ( INOUT ) :: drpot
TYPE ( coeff_type ), TARGET, INTENT ( IN ) :: rhos
INTEGER, DIMENSION ( 3 ), INTENT ( IN ) :: center
REAL ( dbl ), INTENT ( OUT ) :: force
! Locals
INTEGER, DIMENSION ( :, : ), POINTER :: bs
INTEGER, DIMENSION ( : ), POINTER :: nb
INTEGER, DIMENSION ( : ), POINTER :: ns
INTEGER, DIMENSION ( 3 ) :: nc
INTEGER :: i, j, k, ii, ia
LOGICAL :: folded
!------------------------------------------------------------------------------
bs => rhos % pw % pw_grid % bounds
nb => drpot % npts
ns => rhos % pw % pw_grid % npts
folded = .FALSE.
DO i = bs ( 1, 1 ), bs ( 2, 1 )
ia = i - bs ( 1, 1 ) + 1
ii = center ( 1 ) + i - drpot % lb_local ( 1 )
IF ( ii < 0 ) THEN
drpot % px ( ia ) = ii + nb ( 1 ) + 1
folded = .TRUE.
ELSEIF ( ii >= nb ( 1 ) ) THEN
drpot % px ( ia ) = ii - nb ( 1 ) + 1
folded = .TRUE.
ELSE
drpot % px ( ia ) = ii + 1
ENDIF
END DO
DO i = bs ( 1, 2 ), bs ( 2, 2 )
ia = i - bs ( 1, 2 ) + 1
ii = center ( 2 ) + i - drpot % lb_local ( 2 )
IF ( ii < 0 ) THEN
drpot % py ( ia ) = ii + nb ( 2 ) + 1
folded = .TRUE.
ELSEIF ( ii >= nb ( 2 ) ) THEN
drpot % py ( ia ) = ii - nb ( 2 ) + 1
folded = .TRUE.
ELSE
drpot % py ( ia ) = ii + 1
ENDIF
END DO
DO i = bs ( 1, 3 ), bs ( 2, 3 )
ia = i - bs ( 1, 3 ) + 1
ii = center ( 3 ) + i - drpot % lb_local ( 3 )
IF ( ii < 0 ) THEN
drpot % pz ( ia ) = ii + nb ( 3 ) + 1
folded = .TRUE.
ELSEIF ( ii >= nb ( 3 ) ) THEN
drpot % pz ( ia ) = ii - nb ( 3 ) + 1
folded = .TRUE.
ELSE
drpot % pz ( ia ) = ii + 1
ENDIF
END DO
IF ( folded ) THEN
CALL dg_int_patch_1d ( drpot % r, rhos % pw % cr3d, force, ns, &
drpot % px, drpot % py, drpot % pz )
ELSE
nc ( 1 ) = drpot % px ( 1 ) - 1
nc ( 2 ) = drpot % py ( 1 ) - 1
nc ( 3 ) = drpot % pz ( 1 ) - 1
CALL dg_int_patch_1d ( drpot % r, rhos % pw % cr3d, force, ns, nc )
END IF
END SUBROUTINE dg_sum_patch_force_coef_1d
!******************************************************************************
@ -740,7 +1011,7 @@ SUBROUTINE dg_get_patch_1 ( rho0, rhos1, charge1, ex1, ey1, ez1 )
n1 = SIZE ( cd, 1 )
n2 = SIZE ( cd, 2 )
n3 = SIZE ( cd, 3 )
IF ( n1 /= nd(1) .OR. n1 /= nd(1) .OR. n1 /= nd(1) ) THEN
IF ( n1 /= nd ( 1 ) .OR. n2 /= nd ( 2 ) .OR. n3 /= nd ( 3 ) ) THEN
DEALLOCATE ( zs, STAT = ierr )
IF ( ierr /= 0 ) call stop_memory ( "dg_get_patch","zs")
DEALLOCATE ( cd, STAT = ierr )
@ -748,20 +1019,23 @@ SUBROUTINE dg_get_patch_1 ( rho0, rhos1, charge1, ex1, ey1, ez1 )
ENDIF
ENDIF
IF ( .NOT. ALLOCATED ( cd ) ) THEN
ALLOCATE ( zs ( nd(1)*nd(2) ), STAT = ierr )
ALLOCATE ( zs ( nd ( 1 ) * nd ( 2 ) ), STAT = ierr )
IF ( ierr /= 0 ) call stop_memory ( "dg_get_patch","zs",nd(1)*nd(2))
zs = 0._dbl
ALLOCATE ( cd ( nd(1),nd(2),nd(3) ), STAT = ierr )
IF ( ierr /= 0 ) call stop_memory ( "dg_get_patch","cd",n)
cd = 0._dbl
zs=0._dbl
ALLOCATE ( cd ( nd ( 1 ), nd ( 2 ), nd ( 3 ) ), STAT = ierr )
IF ( ierr /= 0 ) call stop_memory ( "dg_get_patch", "cd", n )
cd=0._dbl
ENDIF
za = CMPLX ( 0._dbl, 0._dbl, dbl )
zb = CMPLX ( charge1, 0._dbl, dbl )
CALL rankup ( nd, za, cd, zb, ex1, ey1, ez1, zs )
CALL vr_x_vc ( n, rho0 % pw % cr3d, cd )
IF ( rho0 % pw % in_use == PW_REALDATA3D ) &
CALL vr_x_vc ( n, rho0 % pw % cr3d, cd )
IF ( rho0 % pw % in_use == PW_COMPLEXDATA3D ) &
CALL vc_x_vc ( n, rho0 % pw % cc3d, cd )
CALL fft3d ( BWFFT, nd, cd )
CALL copy_cr ( n, cd, rhos1%pw%cr3d )
CALL copy_cr ( n, cd, rhos1 % pw % cr3d )
END SUBROUTINE dg_get_patch_1
@ -795,7 +1069,7 @@ SUBROUTINE dg_get_patch_2 ( rho0, rhos1, rhos2, charge1, charge2, &
n1 = SIZE ( cd, 1 )
n2 = SIZE ( cd, 2 )
n3 = SIZE ( cd, 3 )
IF ( n1 /= nd(1) .OR. n1 /= nd(1) .OR. n1 /= nd(1) ) THEN
IF ( n1 /= nd ( 1 ) .OR. n2 /= nd ( 2 ) .OR. n3 /= nd ( 3 ) ) THEN
DEALLOCATE ( zs, STAT = ierr )
IF ( ierr /= 0 ) call stop_memory ( "dg_get_patch","zs")
DEALLOCATE ( cd, STAT = ierr )
@ -805,10 +1079,10 @@ SUBROUTINE dg_get_patch_2 ( rho0, rhos1, rhos2, charge1, charge2, &
IF ( .NOT. ALLOCATED ( cd ) ) THEN
ALLOCATE ( zs ( nd(1)*nd(2) ), STAT = ierr )
IF ( ierr /= 0 ) call stop_memory ( "dg_get_patch","zs",nd(1)*nd(2))
zs = 0._dbl
zs=0._dbl
ALLOCATE ( cd ( nd(1),nd(2),nd(3) ), STAT = ierr )
IF ( ierr /= 0 ) call stop_memory ( "dg_get_patch","cd",n)
cd = 0._dbl
cd=0._dbl
ENDIF
za = CMPLX ( 0._dbl, 0._dbl, dbl )
@ -817,9 +1091,12 @@ SUBROUTINE dg_get_patch_2 ( rho0, rhos1, rhos2, charge1, charge2, &
za = CMPLX ( 0._dbl, 1._dbl, dbl )
zb = CMPLX ( charge1, 0._dbl, dbl )
CALL rankup ( nd, za, cd, zb, ex1, ey1, ez1, zs )
CALL vr_x_vc ( n, rho0 % pw % cr3d, cd )
IF ( rho0 % pw % in_use == PW_REALDATA3D ) &
CALL vr_x_vc ( n, rho0 % pw % cr3d, cd )
IF ( rho0 % pw % in_use == PW_COMPLEXDATA3D ) &
CALL vc_x_vc ( n, rho0 % pw % cc3d, cd )
CALL fft3d ( BWFFT, nd, cd )
CALL copy_cri ( n, cd, rhos1%pw%cr3d, rhos2%pw%cr3d )
CALL copy_cri ( n, cd, rhos1 % pw % cr3d, rhos2 % pw % cr3d )
END SUBROUTINE dg_get_patch_2
@ -884,7 +1161,7 @@ END SUBROUTINE dg_add_patch_folded
!******************************************************************************
SUBROUTINE dg_int_patch_simple ( rbx, rby, rbz, rs, f, ns, nc )
SUBROUTINE dg_int_patch_simple_3d ( rbx, rby, rbz, rs, f, ns, nc )
IMPLICIT NONE
@ -914,11 +1191,44 @@ SUBROUTINE dg_int_patch_simple ( rbx, rby, rbz, rs, f, ns, nc )
END DO
END DO
END SUBROUTINE dg_int_patch_simple
END SUBROUTINE dg_int_patch_simple_3d
!******************************************************************************
SUBROUTINE dg_int_patch_folded ( rbx, rby, rbz, rs, f, ns, px, py, pz )
SUBROUTINE dg_int_patch_simple_1d ( rb, rs, f, ns, nc )
IMPLICIT NONE
! Arguments
REAL ( dbl ), DIMENSION (:,:,:), INTENT ( IN ) :: rb
REAL ( dbl ), DIMENSION (:,:,:), INTENT ( IN ) :: rs
REAL ( dbl ), INTENT ( OUT ) :: f
INTEGER, INTENT ( IN ) :: ns ( 3 ), nc ( 3 )
! LOCAL
INTEGER :: i, j, k, ii, jj, kk
REAL ( dbl ) :: s
!------------------------------------------------------------------------------
f = 0._dbl
DO k = 1, ns ( 3 )
kk = nc ( 3 ) + k
DO j = 1, ns ( 2 )
jj = nc ( 2 ) + j
DO i = 1, ns ( 1 )
ii = nc ( 1 ) + i
s = rs ( i, j, k )
f = f + s * rb ( ii, jj, kk )
END DO
END DO
END DO
END SUBROUTINE dg_int_patch_simple_1d
!******************************************************************************
SUBROUTINE dg_int_patch_folded_3d ( rbx, rby, rbz, rs, f, ns, px, py, pz )
IMPLICIT NONE
@ -949,7 +1259,41 @@ SUBROUTINE dg_int_patch_folded ( rbx, rby, rbz, rs, f, ns, px, py, pz )
END DO
END DO
END SUBROUTINE dg_int_patch_folded
END SUBROUTINE dg_int_patch_folded_3d
!******************************************************************************
SUBROUTINE dg_int_patch_folded_1d ( rb, rs, f, ns, px, py, pz )
IMPLICIT NONE
! Arguments
REAL ( dbl ), DIMENSION (:,:,:), INTENT ( IN ) :: rb
REAL ( dbl ), DIMENSION (:,:,:), INTENT ( IN ) :: rs
REAL ( dbl ), INTENT ( INOUT ) :: f
INTEGER, INTENT ( IN ) :: ns ( 3 )
INTEGER, DIMENSION ( : ), INTENT ( IN ) :: px, py, pz
! LOCAL
INTEGER :: i, j, k, ii, jj, kk
REAL ( dbl ) :: s
!------------------------------------------------------------------------------
f = 0._dbl
DO k = 1, ns ( 3 )
kk = pz ( k )
DO j = 1, ns ( 2 )
jj = py ( j )
DO i = 1, ns ( 1 )
ii = px ( i )
s = rs ( i, j, k )
f = f + s * rb ( ii, jj, kk )
END DO
END DO
END DO
END SUBROUTINE dg_int_patch_folded_1d
!******************************************************************************

10
src/ewalds.F
View file

@ -21,7 +21,7 @@
MODULE ewalds
USE coefficient_types, ONLY : coeff_type
USE dgs, ONLY : dg_grid_setup
USE dgs, ONLY : dg_pme_grid_setup
USE dg_types, ONLY : dg_type
USE dg_rho0s, ONLY : dg_rho0_setup
USE ewald_parameters_types, ONLY : ewald_parameters_type
@ -32,7 +32,7 @@ MODULE ewalds
USE molecule_types, ONLY : particle_node_type
USE message_passing, ONLY : mp_sum
USE particle_types, ONLY : particle_type
USE particle_lists, ONLY : particle_list_select
USE particle_lists, ONLY : particle_list_ewald
USE pw_grid_types, ONLY : pw_grid_type, HALFSPACE
USE pw_grids, ONLY : pw_find_cutoff, pw_grid_setup, pw_grid_construct
USE realspace_grid_types, ONLY : realspace_grid_type, rs_grid_setup
@ -518,7 +518,7 @@ SUBROUTINE ewald_initialize ( dg, part, pnode, pnode_grp, ewald_param, box, &
npts_s ( : ) = ewald_param % ns_max
CALL dg_grid_setup ( box, npts_s, ewald_param % epsilon, &
CALL dg_pme_grid_setup ( box, npts_s, ewald_param % epsilon, &
ewald_param % alpha, pme_small_grid, &
pme_big_grid, ewald_param % ewald_type, pnode_grp )
@ -526,7 +526,7 @@ SUBROUTINE ewald_initialize ( dg, part, pnode, pnode_grp, ewald_param, box, &
pnode_grp % scr )
! set up a list with all particles having a chage
CALL particle_list_select ( pnode, dg % plist, dg % nparts )
CALL particle_list_ewald ( pnode, dg % plist, dg % nparts )
IF ( pnode_grp % ionode ) THEN
WRITE ( iw, '( A,T71,E10.4 )' ) &
@ -556,7 +556,7 @@ SUBROUTINE ewald_initialize ( dg, part, pnode, pnode_grp, ewald_param, box, &
CALL rs_grid_setup ( rs, spme_grid, npts_s, pnode_grp % scr )
! set up a list with all particles having a chage
CALL particle_list_select ( pnode, dg % plist, dg % nparts )
CALL particle_list_ewald ( pnode, dg % plist, dg % nparts )
END IF

44
src/extended_system_types.F Normal file
View file

@ -0,0 +1,44 @@
!-----------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!***** cp2k/extended_system_types [1.0] *
!!
!! NAME
!! extended_system_types
!!
!! FUNCTION
!! Lumps all possible extended system variables into one
!! type for easy access and passing
!!
!! AUTHOR
!! CJM
!!
!! MODIFICATION HISTORY
!!
!! SOURCE
!******************************************************************************
MODULE extended_system_types
USE kinds, ONLY : dbl
USE nose, ONLY : lnhc_parameters_type, npt_info_type
IMPLICIT NONE
PRIVATE
PUBLIC :: extended_system_type
! PUBLIC :: npt_info, nhc_part, nhc_baro, nhc_pol
TYPE extended_system_type
TYPE ( lnhc_parameters_type ), POINTER :: nhc_part ( : )
TYPE ( lnhc_parameters_type ), POINTER :: nhc_baro ( : )
TYPE ( lnhc_parameters_type ) :: nhc_pol
TYPE ( npt_info_type ), DIMENSION ( :, : ), POINTER :: npt_info
END TYPE extended_system_type
END MODULE extended_system_types
!!*****
!******************************************************************************

1
src/fftw_lib.F
View file

@ -235,7 +235,6 @@ SUBROUTINE fft3d ( fsign, scale, n, zg, zg_out )
n1 = n(1)
n2 = n(2)
n3 = n(3)
IF ( PRESENT ( zg_out ) ) THEN
fft_in_place = .false.
zgout => zg_out

1
src/fist_nonbond_force.F
View file

@ -312,7 +312,6 @@ SUBROUTINE bonded_correct_gaussian ( ewald_param, molecule, v_bonded_corr, &
v_bonded_corr = 0.0_dbl
!
! return !wdbg
!
alpha = ewald_param % alpha
const = 2.0_dbl * alpha / SQRT( pi )

86
src/fit_force_control.F Normal file
View file

@ -0,0 +1,86 @@
!-----------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/force_control [1.0] *
!!
!! NAME
!! fit_force_control
!!
!! FUNCTION
!! Interface routine for the fit calculations
!!
!! SOURCE
!******************************************************************************
MODULE fit_force_control
USE ewald_parameters_types, ONLY : ewald_parameters_type
USE fist_force, ONLY : fist_force_control => force_control
USE global_types, ONLY : global_environment_type
USE kinds, ONLY : dbl
USE md, ONLY : simulation_parameters_type, thermodynamic_type
USE structure_types, ONLY : structure_type
USE pair_potential, ONLY : potentialparm_type
USE fit_parameters_types, ONLY : fit_parameter_type
USE termination, ONLY : stop_program
IMPLICIT NONE
PRIVATE
PUBLIC :: fit_force
!!*****
!******************************************************************************
CONTAINS
!******************************************************************************
!!****** force_control/force [1.0] *
!!
!! NAME
!! force
!!
!******************************************************************************
SUBROUTINE fit_force ( struc, potparm, thermo, ewald_param, &
fit_param, globenv)
IMPLICIT NONE
! Arguments
TYPE ( structure_type ), INTENT ( INOUT ) :: struc
TYPE ( thermodynamic_type ), INTENT ( INOUT ) :: thermo
TYPE (potentialparm_type), DIMENSION (:,:), INTENT (INOUT) :: potparm
TYPE ( ewald_parameters_type ), INTENT ( INOUT ) :: ewald_param
TYPE (fit_parameter_type), INTENT(INOUT) :: fit_param
TYPE ( global_environment_type), INTENT ( IN ) :: globenv
!------------------------------------------------------------------------------
!locals
LOGICAL, PARAMETER :: box_change = .FALSE.
thermo % ptens_aval = .FALSE.
SELECT CASE ( fit_param % fit_type )
CASE DEFAULT
CALL stop_program ( 'fit_force', 'not implemented' )
CASE ( 'AMOEBA_RI' )
CALL fist_force_control ( struc % molecule, struc % pnode, struc % part, &
struc % box, struc % box_ref, thermo, potparm, ewald_param, &
box_change, globenv, struc%ll_data (1) )
CASE ( 'AMOEBA_POL' )
CALL stop_program ( 'fit_force', 'not implemented' )
END SELECT
END SUBROUTINE fit_force
!!*****
!******************************************************************************
END MODULE fit_force_control
!******************************************************************************

190
src/fit_input.F Normal file
View file

@ -0,0 +1,190 @@
!-----------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 1999 MPI fuer Festkoerperforschung, Stuttgart !
!-----------------------------------------------------------------------------!
MODULE fit_input
USE global_types, ONLY : global_environment_type
USE kinds, ONLY : dbl
USE fit_parameters_types, ONLY : fit_parameter_type
USE parser, ONLY : parser_init, parser_end, read_line, test_next, &
cfield, p_error, get_real, get_int, stop_parser
USE string_utilities, ONLY : uppercase
USE termination, ONLY : stop_program
IMPLICIT NONE
PRIVATE
PUBLIC :: read_fit_section
CONTAINS
!******************************************************************************
!!>---------------------------------------------------------------------------!
!! SECTION: &fit ... &end !
!! !
!! fit_type: [AMOEBA_RI,AMOEBA_POL,FORCE_POL] type of fit run
!! pot_type: [LJ, WILL] type of two body potential to fit !
!! restart: [INITIAL, ALL] restart option !
!! data type: [POS_FORCE, POS_VEL] type of data to read !
!! out unit: [ATOMIC, KELVIN] output unit !
!! tol: max_tol maximum tolerance (fractional ) !
!! iter: max_iter maximum number of iterations !
!! nsteps: nsteps number of trajectory steps !
!! iskip: iskip analyse every iskip step !
!! idump: idump write to dump file every idump steps !
!! iwrite: iwrite write to screen every iwrite steps !
!! !
!!<---------------------------------------------------------------------------!
SUBROUTINE read_fit_section ( fit_param, globenv )
IMPLICIT NONE
! Arguments
TYPE ( fit_parameter_type ), INTENT ( OUT ) :: fit_param
TYPE ( global_environment_type ), INTENT ( IN ) :: globenv
! Locals
INTEGER :: ierror, ilen, iw
CHARACTER ( len = 20 ) :: string
CHARACTER ( len = 5 ) :: label
!------------------------------------------------------------------------------
iw = globenv % scr
!..defaults
fit_param % max_tol = 1.0E-6_dbl
fit_param % max_iter = 10000
fit_param % nsteps = 1
fit_param % iskip = 1
fit_param % iwrite = 1
fit_param % pot_type = 'LJ'
fit_param % fit_type = 'AMOEBA_RI'
fit_param % datatype = 'POS_VEL'
fit_param % outunit = 'ATOMIC'
fit_param % idump = 10
fit_param % restart = 'INITIAL'
!..parse the input section
IF ( globenv % ionode ) THEN
label = '&FIT'
CALL parser_init ( globenv % input_file_name, label, ierror, globenv )
IF ( ierror /= 0 ) THEN
CALL stop_program ( 'read_fit_section', 'no input section &FIT found' )
END IF
CALL read_line()
DO WHILE ( test_next ( ) /= 'X' )
ilen = 8
CALL cfield ( string, ilen )
CALL uppercase ( string )
SELECT CASE ( string )
CASE DEFAULT
CALL p_error ( )
CALL stop_program ( 'read_fit_section', 'unknown option' )
CASE ( 'FIT_TYPE' )
ilen = 10
CALL cfield ( string, ilen )
CALL uppercase ( string )
fit_param % fit_type = string ( 1:ilen )
CASE ( 'POT_TYPE' )
ilen = 6
CALL cfield ( string, ilen )
CALL uppercase ( string )
fit_param % pot_type = string ( 1:ilen )
CASE ( 'DATA' )
ilen = 10
CALL cfield ( string, ilen )
CALL uppercase ( string )
fit_param % datatype = string ( 1:ilen )
CASE ( 'RESTART' )
ilen = 7
CALL cfield ( string, ilen )
CALL uppercase ( string )
fit_param % restart = string ( 1:ilen )
CASE ( 'OUTUNIT' )
ilen = 6
CALL cfield ( string, ilen )
CALL uppercase ( string )
fit_param % outunit = string ( 1:ilen )
CASE ( 'TOL' )
fit_param % max_tol = get_real ( )
CASE ( 'ITER' )
fit_param % max_iter = get_int ( )
CASE ( 'NSTEPS' )
fit_param % nsteps = get_int ( )
CASE ( 'IDUMP' )
fit_param % idump = get_int ( )
CASE ( 'ISKIP' )
fit_param % iskip = get_int ( )
CASE ( 'IWRITE' )
fit_param % iwrite = get_int ( )
END SELECT
! Get the next line
CALL read_line
END DO
END IF
CALL parser_end
!..end of parsing the input section
!..write some information to output
IF ( globenv % print_level >= 0 ) THEN
WRITE ( iw, '( A )' ) ' FIT| Fit Protocol '
WRITE ( iw, '( A, T71, A )' ) &
' FIT| Fit type ', ADJUSTR ( fit_param % fit_type )
WRITE ( iw, '( A, T71, A )' ) &
' FIT| 2 body potential type ', ADJUSTR ( fit_param % pot_type )
WRITE ( iw, '( A, T71, A )' ) &
' FIT| Data type ', ADJUSTR ( fit_param % datatype )
WRITE ( iw, '( A, T71, A )' ) &
' FIT| Restart type', ADJUSTR ( fit_param % restart )
WRITE ( iw, '( A, T71, A )' ) &
' FIT| Out Units ', ADJUSTR ( fit_param % outunit )
WRITE ( iw, '( A, T71, I10 )' ) &
' FIT| Number of steps ', fit_param % nsteps
WRITE ( iw, '( A, T71, I10 )' ) &
' FIT| Iskip ', fit_param % iskip
WRITE ( iw, '( A, T71, I10 )' ) &
' FIT| Idump ', fit_param % idump
WRITE ( iw, '( A, T71, I10 )' ) &
' FIT| iwrite ', fit_param % iwrite
WRITE ( iw, '( A, T71, G10.4 )' ) &
' FIT| Max tolerance ', fit_param % max_tol
WRITE ( iw, '( A, T71, I10 )' ) &
' FIT| Max iterations ', fit_param % max_iter
END IF
END SUBROUTINE read_fit_section
!******************************************************************************
END MODULE fit_input
!******************************************************************************

473
src/fit_module.F Normal file
View file

@ -0,0 +1,473 @@
!------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 1999 MPI fuer Festkoerperforschung, Stuttgart !
!------------------------------------------------------------------------------!
MODULE fit_module
!-----------------------------------------------------------------------------!
!!****** cp2k/fit_module [1.0] *
!!
!! NAME
!! fit_module
!!
!! FUNCTION
!!
!! enables to do different types of fit (rigid_ion, polarizable ..)
!! using the amoeba procedire
!!
!*****************************************************************
!------------------------------------------------------------------------------!
!
USE amoeba, ONLY : amoeba_evaluate
USE global_types, ONLY : global_environment_type
USE kinds, ONLY : dbl
USE md, ONLY : simulation_parameters_type, thermodynamic_type
USE pair_potential, ONLY : potentialparm_type, spline_nonbond_control
USE parser, ONLY : parser_init, parser_end, read_line, test_next, &
cfield, p_error, get_real, get_int
USE particle_types, ONLY : particle_prop_type
USE structure_types, ONLY : structure_type
USE termination, ONLY : stop_program, stop_memory
USE fit_parameters_types, ONLY : fit_parameter_type, fit_averages_type, &
linklist_parameters_type, fit_info_type
USE fit_utilities, ONLY : force_ref, fit_print_info, get_will_params, &
get_lj_params, fit_read_restart, fit_write_restart
USE fit_types, ONLY : copy_params, initialize_params
USE force_fields, ONLY : ATOMNAMESLENGTH
USE simulation_cell, ONLY : get_hinv, cell_type
USE fit_force_control, ONLY : fit_force
USE ewald_parameters_types, ONLY : ewald_parameters_type
!gt this is to fix to do the polarizable fit
! USE force_control, ONLY : force
! USE basis_kind_info, ONLY : kind_info_type
! USE minimize, ONLY : fit_min_control, dump_coef => dump
! USE pol_lists, ONLY : pol_es_lists
! USE pol_electrostatic_force, ONLY : force_electrostatic
! USE pol_overlap_force, ONLY : force_overlap
! USE pol_classic_force, ONLY : force_classic
IMPLICIT NONE
!
PUBLIC :: control_fit, local
PRIVATE
!
!------------------------------------------------------------------------------!
!
CONTAINS
!
!------------------------------------------------------------------------------!
!
!!>----------------------------------------------------------------------------!
!! SECTION: &fit ... &end !
!! !
!! type: [LJ, WILL] type of fit !
!! restart: [INITIAL, ALL] restart option !
!! data type: [POS_FORCE, POS_VEL] type of data to read !
!! out unit: [ATOMIC, KELVIN] output unit !
!! tol: max_tol maximum tolerance (fractional) !
!! iter: max_iter maximum number of iterations !
!! runtype: [AMOEBA_RI, FORCES_POL
!! nsteps: ntime number of trajectory steps !
!! iskip: iskip analyse every iskip step !
!! idump: idump write to dump file every idump steps !
!! !
!! !
!!<----------------------------------------------------------------------------!
!------------------------------------------------------------------------------!
!
SUBROUTINE control_fit(struc,potparm,simpar,fit_param, pstat, &
ewald_param, atom_names,globenv)
IMPLICIT NONE
TYPE (structure_type), INTENT (INOUT) :: struc
TYPE (potentialparm_type ), INTENT ( INOUT ), DIMENSION ( :, : ) :: potparm
TYPE (simulation_parameters_type), INTENT ( IN ) :: simpar
TYPE (particle_prop_type ), DIMENSION ( : ), POINTER :: pstat
TYPE (fit_parameter_type), INTENT(INOUT) :: fit_param
TYPE ( ewald_parameters_type ), INTENT ( INOUT ) :: ewald_param
TYPE ( global_environment_type ), INTENT ( IN ) :: globenv
CHARACTER ( LEN = ATOMNAMESLENGTH ), DIMENSION ( : ), POINTER :: atom_names
!
INTEGER :: isos
! locals
TYPE ( thermodynamic_type ) :: thermo
TYPE ( fit_averages_type ) :: fit_ave
TYPE ( linklist_parameters_type ), POINTER, DIMENSION ( :, : ) :: param
LOGICAL :: conv_flag
!
SELECT CASE (fit_param%fit_type)
!
! before fitting short-range parameters to the cp forces, correct them
! by subtracting the polarization contribution (and Coulomb if present)
!
CASE ('FORCES_POL')
! gt to be implemanted
! CALL calculate_polar_forces(struc,thermo,simpar,potparm)
! gt to be implemanted
!
! amoeba fit on the cpmd forces , or on the forces obtained
! by subtracting the polarizable forces to the cpmd ones.
! This works both for rigid ions model and polarizable force fields
! (no variation of the coefficients, only amoeba minimizer)
!
CASE ('AMOEBA_RI')
! copy all relevent variables to local, global variables
ALLOCATE (param(size(potparm,1),size(potparm,2)),STAT=isos)
IF (isos/=0) CALL stop_memory('fit_control','param')
! allocating the arrays for the analysis
SELECT CASE (fit_param%pot_type)
CASE ('LJ')
ALLOCATE (fit_ave % ave_sig(size(potparm,1),size(potparm,2)),STAT=isos)
IF ( isos /= 0 ) CALL stop_memory ( 'control_fit', 'ave_sig')
ALLOCATE (fit_ave % ave_eps(size(potparm,1),size(potparm,2)),STAT=isos)
IF (isos/=0) CALL stop_memory('control_fit', 'ave_eps')
CASE ('WILL')
ALLOCATE (fit_ave % ave_a(size(potparm,1),size(potparm,2)),STAT=isos)
IF ( isos /= 0 ) CALL stop_memory ( 'control_fit', 'ave_a')
ALLOCATE (fit_ave % ave_b(size(potparm,1),size(potparm,2)),STAT=isos)
IF ( isos /= 0 ) CALL stop_memory ( 'control_fit', 'ave_b')
ALLOCATE (fit_ave % ave_c(size(potparm,1),size(potparm,2)),STAT=isos)
IF ( isos /= 0 ) CALL stop_memory ( 'control_fit', 'ave_c')
END SELECT
conv_flag = .FALSE.
CALL fit(struc,potparm,thermo,simpar,fit_param,param,pstat,fit_ave, &
ewald_param,atom_names,globenv,conv_flag)
IF (conv_flag) THEN
WRITE (*,*)
WRITE (*,'(T23,A)') '******************************'
WRITE (*,'(T23,A)') '* *'
WRITE (*,'(T23,A)') '* FIT|CONVERGENCE ACHIEVED *'
WRITE (*,'(T23,A)') '* *'
WRITE (*,'(T23,A)') '******************************'
WRITE (*,*)
ELSE
WRITE (*,*)
WRITE (*,'(T23,A)') '******************************'
WRITE (*,'(T23,A)') '* *'
WRITE (*,'(T23,A)') '* FIT|WARNING:NO CONVERGENCE *'
WRITE (*,'(T23,A)') '* *'
WRITE (*,'(T23,A)') '******************************'
WRITE (*,*)
END IF
! deallocating the arrays for the analysis
SELECT CASE (fit_param%pot_type)
CASE ('LJ')
IF ( ASSOCIATED ( fit_ave % ave_sig ) ) THEN
DEALLOCATE (fit_ave % ave_sig , STAT=isos)
IF ( isos /= 0 ) CALL stop_memory ( 'control_fit', 'deallocate ave_sig')
END IF
IF ( ASSOCIATED ( fit_ave % ave_eps ) ) THEN
DEALLOCATE (fit_ave % ave_eps , STAT=isos)
IF ( isos /= 0 ) CALL stop_memory ( 'control_fit', 'deallocate ave_eps')
END IF
CASE ('WILL')
IF ( ASSOCIATED ( fit_ave % ave_a ) ) THEN
DEALLOCATE (fit_ave % ave_a , STAT=isos)
IF ( isos /= 0 ) CALL stop_memory ( 'control_fit', 'deallocate ave_a' )
END IF
IF ( ASSOCIATED ( fit_ave % ave_b ) ) THEN
DEALLOCATE (fit_ave % ave_b , STAT=isos)
IF ( isos /= 0 ) CALL stop_memory ( 'control_fit', 'deallocate ave_b' )
END IF
IF ( ASSOCIATED ( fit_ave % ave_c ) ) THEN
DEALLOCATE (fit_ave % ave_c , STAT=isos)
IF ( isos /= 0 ) CALL stop_memory ( 'control_fit', 'deallocate ave_c' )
END IF
END SELECT
IF ( ASSOCIATED ( param ) ) THEN
DEALLOCATE ( param , STAT=isos)
IF ( isos /= 0 ) CALL stop_memory ( 'control_fit', 'deallocate param' )
END IF
! to perform amoeba fit in the case of polarizable force fields
! in which the amoeba parameters can also depend on the coefficents
! TO BE IMPLEMENTED!!!!
! CASE ('AMOEBA_POL')
END SELECT
END SUBROUTINE control_fit
!
!------------------------------------------------------------------------------!
!
SUBROUTINE fit(struc,potparm,thermo,simpar,fit_param,param,pstat,fit_ave,&
ewald_param,atom_names,globenv,conv_flag)
! Given the reference force, and the derived force, optimize
! paramters with a least squares fit.
IMPLICIT NONE
TYPE (structure_type), INTENT (INOUT) :: struc
TYPE (potentialparm_type), DIMENSION (:,:), INTENT (INOUT) :: potparm
TYPE (thermodynamic_type), INTENT (INOUT) :: thermo
TYPE (simulation_parameters_type), INTENT (IN) :: simpar
TYPE (fit_parameter_type), INTENT(INOUT) :: fit_param
TYPE (linklist_parameters_type ), INTENT(INOUT), DIMENSION ( :, : ) :: param
TYPE ( particle_prop_type ), DIMENSION ( : ), POINTER :: pstat
TYPE (fit_averages_type) :: fit_ave
TYPE ( ewald_parameters_type ), INTENT ( INOUT ) :: ewald_param
TYPE ( global_environment_type ), INTENT ( IN ) :: globenv
CHARACTER ( LEN = ATOMNAMESLENGTH ), DIMENSION ( : ), POINTER :: atom_names
LOGICAL, INTENT (INOUT) :: conv_flag
! locals
INTEGER :: i, j, iter
REAL (dbl) :: tol
REAL ( dbl ), DIMENSION ( : ), POINTER :: y
REAL ( dbl ), DIMENSION ( :, : ), POINTER :: p
TYPE (fit_info_type), DIMENSION (:,:), POINTER :: finfo
TYPE (cell_type) , POINTER, DIMENSION (:) :: box_ref
INTERFACE
FUNCTION func(x,struc,potparm,thermo,simpar,fit_param,box_ref, &
finfo,param,pstat,ewald_param, globenv)
USE kinds, ONLY : dbl
USE structure_types, ONLY : structure_type
USE pair_potential, ONLY : potentialparm_type
USE particle_types, ONLY : particle_prop_type
USE md, ONLY : simulation_parameters_type, thermodynamic_type
USE simulation_cell, ONLY : cell_type
USE fit_parameters_types, ONLY : fit_parameter_type, linklist_parameters_type, &
fit_info_type
USE ewald_parameters_types, ONLY : ewald_parameters_type
USE global_types, ONLY : global_environment_type
IMPLICIT NONE
TYPE (structure_type), INTENT (INOUT) :: struc
TYPE (potentialparm_type), DIMENSION (:,:), INTENT (INOUT) :: potparm
TYPE (thermodynamic_type), INTENT (INOUT) :: thermo
TYPE (simulation_parameters_type), INTENT (IN) :: simpar
REAL (dbl), INTENT (IN), DIMENSION (:) :: x
TYPE (fit_parameter_type), INTENT(IN) :: fit_param
TYPE (cell_type) , POINTER, DIMENSION (:) :: box_ref
TYPE (fit_info_type), DIMENSION (:,:), POINTER :: finfo
TYPE (linklist_parameters_type ), INTENT(INOUT), DIMENSION ( :, : ) :: param
TYPE (particle_prop_type ), DIMENSION ( : ), POINTER :: pstat
TYPE ( ewald_parameters_type ), INTENT ( INOUT ) :: ewald_param
TYPE ( global_environment_type), INTENT ( IN ) :: globenv
REAL (dbl) :: func
END FUNCTION func
END INTERFACE
! Allocate and intialize arrays for amoeba
SELECT CASE (fit_param%fit_type)
CASE ('AMOEBA_RI')
CALL initialize_params(fit_param, potparm, param, y, p )
CASE ('AMOEBA_POL')
CALL stop_program("fit", "not implemanted")
END SELECT
! Get reference force
CALL force_ref(struc,potparm,simpar,fit_param,finfo,box_ref)
! initialize y
SELECT CASE ( fit_param % restart )
CASE ( 'ALL' )
CALL fit_read_restart ( y, p )
CASE ( 'INITIAL' )
DO i = 1, SIZE ( y )
y ( i ) = func ( p ( i, : ), struc, potparm, thermo, simpar, &
fit_param, box_ref, finfo, param, pstat, ewald_param, globenv )
END DO
CASE DEFAULT
CALL stop_program ( "fit", "no valid restart option" )
END SELECT
! optimize
DO iter = 1, fit_param%max_iter, fit_param%iwrite
CALL amoeba_evaluate(p,y,tol,struc, thermo, potparm, simpar, fit_param, &
box_ref, finfo, param, pstat, ewald_param, globenv, func )
SELECT CASE (fit_param%fit_type)
CASE ('AMOEBA_RI')
SELECT CASE (fit_param%pot_type)
CASE ('LJ')
CALL get_lj_params(potparm, y, param, fit_ave%ave_sig, fit_ave%ave_eps)
CASE ('WILL')
CALL get_will_params(potparm, y, param, fit_ave% ave_a, &
fit_ave%ave_b, fit_ave%ave_c)
END SELECT
CALL fit_print_info(tol,potparm,fit_param,y,atom_names, &
fit_ave,globenv,iter)
CASE ('AMOEBA_POL')
END SELECT
! write restart file
IF ( MOD ( iter, fit_param % idump ) == 0 ) &
CALL fit_write_restart ( y, p )
! convergence acheived
IF (tol<=fit_param%max_tol) THEN
conv_flag = .TRUE.
CALL fit_write_restart ( y, p )
EXIT
END IF
END DO
END SUBROUTINE fit
!
!------------------------------------------------------------------------------!
!
SUBROUTINE forcefit(struc, potparm, thermo,box_ref,finfo, &
fit_param, ewald_param, globenv)
IMPLICIT NONE
TYPE (structure_type), INTENT (INOUT) :: struc
TYPE (potentialparm_type), DIMENSION (:,:), INTENT (INOUT) :: potparm
TYPE (thermodynamic_type), INTENT (INOUT) :: thermo
TYPE (cell_type) , POINTER, DIMENSION (:) :: box_ref
TYPE (fit_info_type), DIMENSION (:,:), POINTER :: finfo
TYPE (fit_parameter_type), INTENT(INOUT) :: fit_param
TYPE ( ewald_parameters_type ), INTENT ( INOUT ) :: ewald_param
TYPE ( global_environment_type), INTENT ( IN ) :: globenv
!locals
INTEGER :: itime, ipart, npart
! From the reference trajectory file obtain the fit force
! using a CP2K force field.
npart = size(struc%pnode)
DO itime = 1, fit_param%nsteps
DO ipart = 1, npart
struc%pnode(ipart) %p%r(:) = finfo(itime,ipart) %r(:)
END DO
!pass the box
struc%box=box_ref(itime)
CALL fit_force ( struc, potparm, thermo, &
ewald_param, fit_param, globenv )
DO ipart = 1, npart
finfo(itime,ipart) %f_fit(:) = struc%pnode(ipart) %p%f(:)
END DO
END DO
END SUBROUTINE forcefit
!
!------------------------------------------------------------------------------!
!
FUNCTION local(x,struc,potparm,thermo,simpar,fit_param,box_ref,finfo, &
param, pstat, ewald_param, globenv )
IMPLICIT NONE
REAL (dbl), DIMENSION (:), INTENT (IN) :: x
TYPE (structure_type), INTENT (INOUT) :: struc
TYPE (thermodynamic_type), INTENT (INOUT) :: thermo
TYPE (simulation_parameters_type), INTENT (IN) :: simpar
TYPE (potentialparm_type), DIMENSION (:,:), INTENT (INOUT) :: potparm
TYPE (fit_parameter_type), INTENT(INOUT) :: fit_param
TYPE (cell_type) , POINTER, DIMENSION (:) :: box_ref
TYPE (fit_info_type), DIMENSION (:,:), POINTER :: finfo
TYPE (linklist_parameters_type ), INTENT(INOUT), DIMENSION ( :, : ) :: param
TYPE (particle_prop_type ), DIMENSION ( : ), POINTER :: pstat
TYPE ( ewald_parameters_type ), INTENT ( INOUT ) :: ewald_param
TYPE ( global_environment_type), INTENT ( IN ) :: globenv
! locals
REAL (dbl), DIMENSION (3) :: vec
REAL (dbl) :: local
INTEGER :: itime, ipart, npart
npart = size(struc%pnode)
SELECT CASE (fit_param%fit_type)
CASE ('AMOEBA_RI')
CALL copy_params(fit_param, potparm , pstat, param, ewald_param, x)
CALL forcefit(struc,potparm,thermo,box_ref,finfo, &
fit_param, ewald_param, globenv)
CASE default
CALL stop_program('local', 'not implemented')
END SELECT
local = 0._dbl
DO itime = 1, fit_param%nsteps
DO ipart = 1, npart
vec(:) = finfo(itime,ipart) %f_fit(:) - finfo(itime,ipart) %f_ref(:)
local = local + dot_product(vec,vec)
END DO
END DO
! The addition of the 1 is a safety mechanism so
! that one can define the FRACTIONAL tolerance for
! systems whos func value may be converged at zero.
! It does not effect other converged func values.
! dbg
! local = local/(3._dbl*npart*fitpar%ntime) + 1._dbl
! dbg
local=local/(3._dbl*npart*fit_param%nsteps)
END FUNCTION local
!
!------------------------------------------------------------------------------!
!
END MODULE fit_module
!------------------------------------------------------------------------------!

49
src/fit_parameters_types.F Normal file
View file

@ -0,0 +1,49 @@
!-----------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
MODULE fit_parameters_types
USE kinds, ONLY : dbl
IMPLICIT NONE
PRIVATE
PUBLIC :: linklist_real_type, linklist_parameters_type, &
fit_parameter_type, fit_info_type, fit_averages_type
TYPE linklist_real_type
REAL ( dbl ), POINTER :: real_value
TYPE ( linklist_real_type ), POINTER :: next
END TYPE linklist_real_type
TYPE linklist_parameters_type
TYPE ( linklist_real_type ), POINTER :: lleps
TYPE ( linklist_real_type ), POINTER :: llsig
TYPE ( linklist_real_type ), POINTER :: lla
TYPE ( linklist_real_type ), POINTER :: llb
TYPE ( linklist_real_type ), POINTER :: llc
END TYPE linklist_parameters_type
TYPE fit_parameter_type
CHARACTER ( len = 10 ) :: datatype, outunit, restart
CHARACTER ( len = 10 ) :: fit_type, pot_type
INTEGER :: max_iter, nsteps, iskip, iwrite, idump
REAL ( dbl ) :: max_tol
END TYPE fit_parameter_type
TYPE fit_info_type
REAL ( dbl ), DIMENSION ( 3 ) :: f_ref, f_fit, r
END TYPE fit_info_type
TYPE fit_averages_type
REAL (dbl), POINTER, DIMENSION(:,:) :: ave_sig, ave_eps
REAL (dbl), POINTER, DIMENSION(:,:) :: ave_a, ave_b, ave_c
END TYPE fit_averages_type
END MODULE fit_parameters_types
!******************************************************************************
!!*****

385
src/fit_types.F Normal file
View file

@ -0,0 +1,385 @@
!------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 1999 MPI fuer Festkoerperforschung, Stuttgart !
!------------------------------------------------------------------------------!
MODULE fit_types
!-----------------------------------------------------------------------------!
!!****** cp2k/fit_types [1.0] *
!!
!! NAME
!! fit_types
!!
!! FUNCTION
!!
!! prepares the array for the amoeba for different types of potentials
!!
!*****************************************************************
!------------------------------------------------------------------------------!
!
USE amoeba, ONLY : amoeba_evaluate
USE global_types, ONLY : global_environment_type
USE kinds, ONLY : dbl
USE md, ONLY : simulation_parameters_type, thermodynamic_type
USE ewald_parameters_types, ONLY : ewald_parameters_type
USE pair_potential, ONLY : potentialparm_type, spline_nonbond_control
USE parser, ONLY : parser_init, parser_end, read_line, test_next, &
cfield, p_error, get_real, get_int
USE particle_types, ONLY : particle_prop_type
USE structure_types, ONLY : structure_type
USE termination, ONLY : stop_program, stop_memory
USE fit_parameters_types, ONLY : fit_parameter_type, linklist_real_type, &
linklist_parameters_type
IMPLICIT NONE
!
PUBLIC :: copy_params, initialize_params
PRIVATE
!------------------------------------------------------------------------------!
!
CONTAINS
!
!------------------------------------------------------------------------------!
!
SUBROUTINE copy_params(fit_param, potparm, pstat, param, ewald_param, x)
IMPLICIT NONE
TYPE ( fit_parameter_type ), INTENT ( IN ) :: fit_param
TYPE (potentialparm_type), DIMENSION (:,:), INTENT (INOUT) :: potparm
TYPE (linklist_parameters_type ), INTENT(INOUT), DIMENSION ( :, : ) :: param
TYPE ( ewald_parameters_type ), INTENT ( IN ) :: ewald_param
TYPE ( particle_prop_type ), DIMENSION ( : ), POINTER :: pstat
REAL (dbl), DIMENSION (:), INTENT (IN) :: x
!------------------------------------------------------------------------------!
SELECT CASE (fit_param%pot_type)
CASE ('LJ')
CALL lj_fit(potparm,pstat,param,ewald_param,x)
CASE ('WILL')
CALL will_fit(potparm,pstat,param,ewald_param,x)
CASE DEFAULT
CALL stop_program ( "copy_params", "no valid potential type" )
END SELECT
END SUBROUTINE copy_params
!
!------------------------------------------------------------------------------!
!
SUBROUTINE initialize_params(fit_param, potparm, param, y, p)
IMPLICIT NONE
TYPE ( fit_parameter_type ), INTENT ( IN ) :: fit_param
TYPE (potentialparm_type), DIMENSION (:,:), INTENT (INOUT) :: potparm
TYPE (linklist_parameters_type ), INTENT(INOUT), DIMENSION ( :, : ) :: param
REAL ( dbl ), DIMENSION ( : ), POINTER :: y
REAL ( dbl ), DIMENSION ( :, : ), POINTER :: p
!------------------------------------------------------------------------------!
SELECT CASE (fit_param%pot_type)
CASE ('LJ')
CALL lj_initialize(potparm,param,y,p)
CASE ('WILL')
CALL will_initialize(potparm,param,y,p)
CASE DEFAULT
CALL stop_program ( "copy_params", "no valid potential type" )
END SELECT
END SUBROUTINE initialize_params
!
!------------------------------------------------------------------------------!
!
SUBROUTINE lj_fit(potparm,pstat,param,ewald_param,x)
IMPLICIT NONE
TYPE (potentialparm_type), DIMENSION (:,:), INTENT (INOUT) :: potparm
TYPE ( particle_prop_type ), DIMENSION ( : ), POINTER :: pstat
TYPE (linklist_parameters_type ), INTENT(INOUT), DIMENSION ( :, : ) :: param
TYPE ( ewald_parameters_type ), INTENT ( IN ) :: ewald_param
REAL (dbl), DIMENSION (:), INTENT (IN) :: x
! locals
INTEGER :: i, icount, iatom, jatom, numpar, isos
icount = 0
DO iatom = 1, size(potparm,1)
DO jatom = iatom, size(potparm,2)
icount = icount + 1
potparm(iatom,jatom) %lj%epsilon = x(icount)*x(icount)
! symmetrizing
potparm(jatom,iatom) %lj%epsilon = potparm(iatom,jatom) &
%lj%epsilon
icount = icount + 1
potparm(iatom,jatom) %lj%sigma6 = x(icount)*x(icount)
potparm(iatom,jatom) %lj%sigma12 = x(icount)*x(icount) &
*x(icount)*x(icount)
! symmetrizing
potparm(jatom,iatom) %lj%sigma6 = potparm(iatom,jatom) &
%lj%sigma6
potparm(jatom,iatom) %lj%sigma12 = potparm(iatom,jatom) &
%lj%sigma12
END DO
END DO
! initialize splines
CALL spline_nonbond_control(potparm, pstat, 2500, ewald_param)
END SUBROUTINE lj_fit
!
!------------------------------------------------------------------------------!
!
SUBROUTINE will_fit(potparm,pstat,param,ewald_param,x)
IMPLICIT NONE
TYPE (potentialparm_type), DIMENSION (:,:), INTENT (INOUT) :: potparm
TYPE ( particle_prop_type ), DIMENSION ( : ), POINTER :: pstat
TYPE (linklist_parameters_type ), INTENT(INOUT), DIMENSION ( :, : ) :: param
TYPE ( ewald_parameters_type ), INTENT ( IN ) :: ewald_param
REAL (dbl), DIMENSION (:), INTENT (IN) :: x
! locals
INTEGER :: i, icount, iatom, jatom, numpar, isos
icount = 0
DO iatom = 1, size(potparm,1)
DO jatom = iatom, size(potparm,2)
icount = icount + 1
potparm(iatom,jatom) %willis%a = x(icount)*x(icount)
! symmetrizing
potparm(jatom,iatom) %willis%a = potparm(iatom,jatom) &
%willis%a
icount = icount + 1
potparm(iatom,jatom) %willis%b = x(icount)*x(icount)
! symmetrizing
potparm(jatom,iatom) %willis%b = potparm(iatom,jatom) &
%willis%b
icount = icount + 1
potparm(iatom,jatom) %willis%c = x(icount)*x(icount)
! symmetrizing
potparm(jatom,iatom) %willis%c = potparm(iatom,jatom) &
%willis%c
END DO
END DO
! initialize splines
CALL spline_nonbond_control(potparm,pstat,2500,ewald_param)
END SUBROUTINE will_fit
!
!------------------------------------------------------------------------------!
!
SUBROUTINE lj_initialize(potparm,param,y,p)
IMPLICIT NONE
TYPE (potentialparm_type), DIMENSION (:,:), INTENT (INOUT) :: potparm
TYPE (linklist_parameters_type ), INTENT(INOUT), DIMENSION ( :, : ) :: param
REAL ( dbl ), DIMENSION ( : ), POINTER :: y
REAL ( dbl ), DIMENSION ( :, : ), POINTER :: p
! locals
INTEGER :: i, icount, iatom, jatom, numpar, isos
REAL(dbl) :: sig, eps
TYPE (linklist_real_type), POINTER :: llsig
TYPE (linklist_real_type), POINTER :: lleps
! number of potential parameters to be optimized is 2*number
! of interactions (remember: we only need the upper triangle)
! e.g. for LJ each interaction has 2 parameters, sigma and epsilon
numpar = 2*size(potparm,1)*(size(potparm,1)+1)/2
ALLOCATE (y(numpar+1),STAT=isos)
IF (isos/=0) CALL stop_memory('lj_initialize','alloc y')
ALLOCATE (p(numpar+1,numpar),STAT=isos)
IF (isos/=0) CALL stop_memory('lj_initialize','alloc p')
! i=0,numpar is the numpar+1 vector in p(i,j). It is also
! used to generate the initial guesses for the verticies
! of the amoeba
ALLOCATE (lleps,STAT=isos)
IF (isos/=0) CALL stop_memory('lj_initialize', 'lleps')
ALLOCATE (llsig,STAT=isos)
IF (isos/=0) CALL stop_memory('lj_initialize', 'llsig')
icount = 0
DO iatom = 1, size(potparm,1)
DO jatom = iatom, size(potparm,2)
param(iatom,jatom) %lleps => lleps
param(iatom,jatom) %llsig => llsig
icount = icount + 1
! loop for eps
DO i = 1, numpar + 1
eps=potparm(iatom,jatom) %lj%epsilon
p(i,icount) = SQRT(eps)
lleps%real_value => p(i,icount)
ALLOCATE (lleps%next,STAT=isos)
IF (isos/=0) CALL stop_memory('lj_initialize', 'lleps%next' )
lleps => lleps%next
END DO
icount = icount + 1
! loop for sig
DO i = 1, numpar + 1
sig=potparm(iatom,jatom) %lj%sigma6
p(i,icount) = SQRT(sig)
llsig%real_value => p(i,icount)
ALLOCATE (llsig%next,STAT=isos)
IF (isos/=0) CALL stop_memory('lj_initialize', 'llsig%next')
llsig => llsig%next
END DO
END DO
END DO
! scaling the appropriate interaction parameter by 1.1
! to construct the amoeba
DO i = 0, numpar
IF (i>0) p(i+1,i) = p(i+1,i)*1.1_dbl
END DO
END SUBROUTINE lj_initialize
!
!------------------------------------------------------------------------------!
!
SUBROUTINE will_initialize(potparm,param,y,p)
IMPLICIT NONE
TYPE (potentialparm_type), DIMENSION (:,:), INTENT (INOUT) :: potparm
TYPE (linklist_parameters_type ), INTENT(INOUT), DIMENSION ( :, : ) :: param
REAL ( dbl ), DIMENSION ( : ), POINTER :: y
REAL ( dbl ), DIMENSION ( :, : ), POINTER :: p
! locals
INTEGER :: i, icount, iatom, jatom, numpar, isos
REAL(dbl) a,b,c
TYPE (linklist_real_type), POINTER :: lla
TYPE (linklist_real_type), POINTER :: llb
TYPE (linklist_real_type), POINTER :: llc
! of interactions (remember: we only need the upper triangle)
! e.g. for WILL each interaction has 3 parameters, a, b, and c.
numpar = 3_dbl*size(potparm,1)*(size(potparm,1)+1)/2
ALLOCATE (y(numpar+1),STAT=isos)
IF (isos/=0) CALL stop_memory('will_initialize','alloc y')
ALLOCATE (p(numpar+1,numpar),STAT=isos)
IF (isos/=0) CALL stop_memory('will_initialize','alloc p')
! i=0,numpar is the numpar+1 vector in p(i,j). It is also
! used to generate the initial guesses for the verticies
! of the amoeba
ALLOCATE (lla,STAT=isos)
IF (isos/=0) CALL stop_memory('will_initialize','alloc lla')
ALLOCATE (llb,STAT=isos)
IF (isos/=0) CALL stop_memory('will_initialize','alloc llb')
ALLOCATE (llc,STAT=isos)
IF (isos/=0) CALL stop_memory('will_initialize','alloc llc')
icount = 0
DO iatom = 1, size(potparm,1)
DO jatom = iatom, size(potparm,2)
param(iatom,jatom) %lla => lla
param(iatom,jatom) %llb => llb
param(iatom,jatom) %llc => llc
icount = icount + 1
! loop for a
DO i = 1, numpar + 1
a=potparm(iatom,jatom) %willis%a
p(i,icount) = SQRT(a)
lla%real_value => p(i,icount)
ALLOCATE (lla%next,STAT=isos)
IF (isos/=0) CALL stop_memory('will_initialize', 'alloc lla%next')
lla => lla%next
END DO
icount = icount + 1
! loop for b
DO i = 1, numpar + 1
b=potparm(iatom,jatom) %willis%b
p(i,icount) = SQRT(b)
llb%real_value => p(i,icount)
ALLOCATE (llb%next,STAT=isos)
IF (isos/=0) CALL stop_memory('will_initialize', 'alloc llb%next')
llb => llb%next
END DO
icount = icount + 1
! loop for c
DO i = 1, numpar + 1
c=potparm(iatom,jatom) %willis%c
p(i,icount) = SQRT(c)
llc%real_value => p(i,icount)
ALLOCATE (llc%next,STAT=isos)
IF (isos/=0) CALL stop_memory('will_initialize', 'alloc llc%next')
llc => llc%next
END DO
END DO
END DO
! scaling the appropriate interaction parameter by 1.1
! to construct the amoeba
DO i = 0, numpar
IF (i>0) p(i+1,i) = p(i+1,i)*1.1_dbl
END DO
END SUBROUTINE will_initialize
!
!------------------------------------------------------------------------------!
END MODULE fit_types
!-----------------------------------------------------------------------------!

360
src/fit_utilities.F Normal file
View file

@ -0,0 +1,360 @@
!------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 1999 MPI fuer Festkoerperforschung, Stuttgart !
!------------------------------------------------------------------------------!
MODULE fit_utilities
!-----------------------------------------------------------------------------!
!!****** cp2k/fit_utilities [1.0] *
!!
!! NAME
!! fit_utilities
!!
!! FUNCTION
!!
!! enables to do different types of fit (rigid_ion, polarizable ..)
!! using the amoeba procedire
!!
!*****************************************************************
!------------------------------------------------------------------------------!
!
USE global_types, ONLY : global_environment_type
USE kinds, ONLY : dbl
USE md, ONLY : simulation_parameters_type
USE pair_potential, ONLY : potentialparm_type
USE structure_types, ONLY : structure_type
USE termination, ONLY : stop_program, stop_memory
USE fit_parameters_types, ONLY : fit_parameter_type, fit_averages_type, &
linklist_parameters_type, fit_info_type, linklist_real_type
USE force_fields, ONLY : ATOMNAMESLENGTH
USE simulation_cell, ONLY : get_hinv, cell_type
IMPLICIT NONE
!
PUBLIC :: force_ref, fit_print_info, get_will_params, get_lj_params, &
fit_read_restart, fit_write_restart
PRIVATE
!------------------------------------------------------------------------------!
!
CONTAINS
!
!------------------------------------------------------------------------------!
!
SUBROUTINE force_ref(struc,potparm,simpar,fit_param,finfo, box_ref)
IMPLICIT NONE
TYPE (structure_type), INTENT (INOUT) :: struc
TYPE (potentialparm_type), DIMENSION (:,:), INTENT (INOUT) :: potparm
TYPE (simulation_parameters_type), INTENT (IN) :: simpar
TYPE ( fit_parameter_type ), INTENT ( INOUT ) :: fit_param
TYPE ( fit_info_type ), DIMENSION( :, : ), POINTER :: finfo
TYPE (cell_type) , POINTER, DIMENSION (:) :: box_ref
! locals
INTEGER :: npart, ipart, itime, ii, it, isos, jj, kk
REAL (dbl), ALLOCATABLE, DIMENSION (:,:,:) :: r, v
REAL (dbl), DIMENSION (3) :: fref
REAL (dbl), DIMENSION(3,3) :: hmat
LOGICAL :: copy
! From the first principles trajectory file obtain the total force
! on the particle by inverting the velcocity Verlet equations
npart = size(struc%pnode)
itime = fit_param%nsteps/fit_param%iskip
ALLOCATE (box_ref(itime), stat = isos)
IF ( isos /= 0 ) CALL stop_memory ( 'force_ref', 'box_ref')
IF ( fit_param % datatype == 'POS_FORCE' ) THEN
ALLOCATE ( finfo ( itime, npart ), STAT = isos )
IF ( isos /= 0 ) CALL stop_memory ( 'force_ref', 'finfo')
ELSE IF ( fit_param % datatype == 'POS_VEL' ) THEN
! Can only compute force up to nsteps-1 step: thus allocated so
IF ( MOD ( fit_param % nsteps, itime ) == 0 ) THEN
ALLOCATE ( finfo ( itime-1, npart ), STAT = isos )
IF ( isos /= 0 ) CALL stop_memory ( 'force_ref', 'finfo')
ELSE
ALLOCATE ( finfo ( itime, npart ), STAT = isos )
IF ( isos /= 0 ) CALL stop_memory ( 'force_ref', 'finfo')
END IF
END IF
SELECT CASE (fit_param%datatype)
! FOR DATA TYPE FORCES
CASE ('POS_FORCE')
ALLOCATE (r(1,1,3),STAT=isos)
ALLOCATE (v(1,1,3),STAT=isos) ! reads force data
IF (isos/=0) CALL stop_memory('force_ref','failed to allocate r and v')
copy = .FALSE.
it = 0
DO itime = 1, fit_param%nsteps
IF (mod(itime,fit_param%iskip)==0) THEN
copy = .TRUE.
it = it + 1
END IF
! particle position and forces from ab-initio trajectory
DO ipart = 1, npart
READ (250,*) ii, r(1,1,:), v(1,1,:)
IF (copy) THEN
finfo(it,ipart) %r(:) = r(1,1,:)
finfo(it,ipart) %f_ref(:) = v(1,1,:)
END IF
END DO
! read the box
DO ii=1,3
READ (250,*) kk, (hmat(ii,jj),jj=1,3)
ENDDO
IF (copy) THEN
box_ref(it)%hmat=hmat
box_ref(it)%perd=1
CALL get_hinv(box_ref(it))
ENDIF
copy = .FALSE.
END DO
! Reassign fit_param%nsteps=it
fit_param%nsteps = it
DEALLOCATE (r,STAT=isos)
DEALLOCATE (v,STAT=isos)
IF (isos/=0) CALL stop_memory('force_ref', 'r and v')
! FOR DATA TYPE TRAJECTORY
CASE ('POS_VEL')
ALLOCATE (r(2,npart,3),STAT=isos)
ALLOCATE (v(2,npart,3),STAT=isos) ! reads velocity data
IF (isos/=0) CALL stop_memory('force_ref','r and v')
! Read in particles from t=1
DO ipart = 1, npart
READ (200,*) ii, r(1,ipart,:), v(1,ipart,:)
END DO
! initialize counter and logical
copy = .FALSE.
it = 1
DO itime = 2, fit_param%nsteps
! set flag
IF (mod(itime-1,fit_param%iskip)==0) copy = .TRUE.
DO ipart = 1, npart
READ (200,*) ii, r(2,ipart,:), v(2,ipart,:)
! Invert velocity verlet equations
fref(:) = 2._dbl*struc%pnode(ipart) %p%prop%mass* &
(r(2,ipart,:)-r(1,ipart,:)-simpar%dt*v(1,ipart,:))/ &
simpar%dt/simpar%dt
! reassigning r(1) to r(2) and v(1) to v(2) for next force calculation
IF (copy) THEN
finfo(it,ipart) %r(:) = r(1,ipart,:)
finfo(it,ipart) %f_ref(:) = fref(:)
END IF
END DO
! increment the counter
IF ((itime/=fit_param%nsteps) .AND. copy) it = it + 1
! reset logical
copy = .FALSE.
r(1,:,:) = r(2,:,:)
v(1,:,:) = v(2,:,:)
END DO
! Reassign fit_param%nsteps=it-1
fit_param%nsteps = it - 1
DEALLOCATE (r,STAT=isos)
DEALLOCATE (v,STAT=isos)
IF (isos/=0) CALL stop_program('force_ref', &
'failed to deallocate r and v')
END SELECT
END SUBROUTINE force_ref
!
!------------------------------------------------------------------------------!
!
SUBROUTINE fit_print_info(tol,potparm,fit_param,y,atom_names,&
fit_ave,globenv,iter)
IMPLICIT NONE
REAL (dbl), INTENT (IN) :: tol
TYPE (potentialparm_type), DIMENSION (:,:), INTENT (IN) :: potparm
TYPE ( fit_parameter_type ), INTENT ( IN ) :: fit_param
REAL ( dbl ), DIMENSION ( : ), INTENT ( IN ) :: y
CHARACTER ( LEN = ATOMNAMESLENGTH ), DIMENSION ( : ), POINTER :: atom_names
TYPE (fit_averages_type) :: fit_ave
TYPE ( global_environment_type ), INTENT ( IN ) :: globenv
INTEGER :: iter
! local variables
INTEGER :: iatom, jatom, iw
iw = globenv % scr
WRITE (iw,*)
WRITE (iw,'(A,A)') ' **************************************', &
'*****************************************'
WRITE (iw,'(A,T71,A)') ' FIT TYPE = ', &
adjustr(fit_param%fit_type)
WRITE (iw,'(A,T71,A)') ' OUT UNIT = ', &
adjustr(fit_param%outunit)
WRITE (iw,'(A,T71,I10)') ' ITERATION NUMBER = ', iter
WRITE (iw,'(A,T66,E15.7)') ' |(F_REF-F_FIT)|^2 = ', &
sum(y)/size(y)
WRITE (iw,'(A,T66,E15.7)') ' TOLERANCE = ', tol
WRITE (iw,*)
IF (fit_param%pot_type=='LJ') THEN
WRITE (iw,'(T5,A,T15,A,T50,A,T73,A)') 'ATOMI', 'ATOMJ', 'EPSILON', &
'SIGMA'
WRITE (iw,*)
DO iatom = 1, size(potparm,1)
DO jatom = iatom, size(potparm,2)
WRITE (iw,'(T5,A,T15,A,T45,E15.8,T66,E15.8)') atom_names(iatom), &
atom_names(jatom), fit_ave% ave_eps(iatom,jatom), &
fit_ave% ave_sig(iatom,jatom)
END DO
END DO
ELSE IF (fit_param%pot_type=='WILL') THEN
WRITE (iw,'(T5,A,T15,A,T44,A,T60,A,T75,A)') 'ATOMI', 'ATOMJ', 'A', &
'B', 'C'
WRITE (iw,*)
DO iatom = 1, size(potparm,1)
DO jatom = iatom, size(potparm,2)
WRITE (iw,'(T5,A,T15,A,T39,E12.6,T54,E12.6,T69,E12.6)') &
atom_names(iatom), atom_names(jatom), fit_ave%ave_a(iatom,jatom), &
fit_ave% ave_b(iatom,jatom), fit_ave %ave_c(iatom,jatom)
END DO
END DO
END IF
WRITE (iw,'(A,A)') ' **************************************', &
'*****************************************'
END SUBROUTINE fit_print_info
!
!------------------------------------------------------------------------------!
!
SUBROUTINE get_lj_params(potparm, y, param, ave_sig, ave_eps)
IMPLICIT NONE
TYPE (potentialparm_type), DIMENSION (:,:), INTENT (IN) :: potparm
REAL ( dbl ), DIMENSION ( : ), INTENT ( IN ) :: y
TYPE ( linklist_parameters_type ), INTENT(IN), DIMENSION ( :, : ) :: param
REAL ( dbl ), POINTER, DIMENSION (:,:) :: ave_sig
REAL ( dbl ), POINTER, DIMENSION (:,:) :: ave_eps
! locals
INTEGER :: iatom, jatom, i
TYPE (linklist_real_type), POINTER :: lleps
TYPE (linklist_real_type), POINTER :: llsig
ave_sig = 0_dbl
ave_eps = 0_dbl
DO iatom = 1, size(potparm,1)
DO jatom = iatom, size(potparm,2)
lleps => param(iatom,jatom) %lleps
llsig => param(iatom,jatom) %llsig
DO i = 1, size(y)
ave_sig(iatom,jatom) = ave_sig(iatom,jatom) + &
llsig%real_value*llsig%real_value
ave_eps(iatom,jatom) = ave_eps(iatom,jatom) + &
lleps%real_value*lleps%real_value
llsig => llsig%next
lleps => lleps%next
END DO
ave_sig(iatom,jatom) = ave_sig(iatom,jatom)/size(y)
ave_eps(iatom,jatom) = ave_eps(iatom,jatom)/size(y)
END DO
END DO
ave_sig = ave_sig**(1._dbl/6._dbl)
END SUBROUTINE get_lj_params
!
!------------------------------------------------------------------------------!
!
SUBROUTINE get_will_params(potparm, y, param, ave_a, ave_b, ave_c)
IMPLICIT NONE
TYPE (potentialparm_type), DIMENSION (:,:), INTENT (IN) :: potparm
REAL ( dbl ), DIMENSION ( : ), INTENT ( IN ) :: y
TYPE ( linklist_parameters_type ), INTENT(IN), DIMENSION ( :, : ) :: param
REAL ( dbl ), POINTER , DIMENSION (:,:) :: ave_a
REAL ( dbl ), POINTER , DIMENSION (:,:) :: ave_b
REAL ( dbl ), POINTER , DIMENSION (:,:) :: ave_c
! locals
INTEGER :: iatom, jatom, i
TYPE (linklist_real_type), POINTER :: lla
TYPE (linklist_real_type), POINTER :: llb
TYPE (linklist_real_type), POINTER :: llc
ave_a = 0_dbl
ave_b = 0_dbl
ave_c = 0_dbl
DO iatom = 1, size(potparm,1)
DO jatom = iatom, size(potparm,2)
lla => param(iatom,jatom) %lla
llb => param(iatom,jatom) %llb
llc => param(iatom,jatom) %llc
DO i = 1, size(y)
ave_a(iatom,jatom) = ave_a(iatom,jatom) + lla%real_value*lla%real_value
ave_b(iatom,jatom) = ave_b(iatom,jatom) + llb%real_value*llb%real_value
ave_c(iatom,jatom) = ave_c(iatom,jatom) + llc%real_value*llc%real_value
lla => lla%next
llb => llb%next
llc => llc%next
END DO
ave_a(iatom,jatom) = ave_a(iatom,jatom)/size(y)
ave_b(iatom,jatom) = ave_b(iatom,jatom)/size(y)
ave_c(iatom,jatom) = ave_c(iatom,jatom)/size(y)
END DO
END DO
END SUBROUTINE get_will_params
!******************************************************************************
SUBROUTINE fit_write_restart ( y, p )
IMPLICIT NONE
! Arguments
REAL ( dbl ), DIMENSION ( : ), INTENT ( IN ) :: y
REAL ( dbl ), DIMENSION ( :, : ), INTENT ( IN ) :: p
! Locals
INTEGER :: i, j, iochannel = 99
!------------------------------------------------------------------------------
! write to the restart file
OPEN ( UNIT = iochannel, FILE = 'amoeba.rst', STATUS = 'UNKNOWN' )
DO i = 1, SIZE ( y )
WRITE ( iochannel, FMT = * ) y ( i )
END DO
DO i = 1, SIZE ( y )
WRITE ( iochannel, FMT = * ) ( p ( i, j ), j = 1, SIZE ( p, 2 ) )
END DO
CLOSE ( UNIT = iochannel )
END SUBROUTINE fit_write_restart
!******************************************************************************
SUBROUTINE fit_read_restart ( y, p )
IMPLICIT NONE
! Arguments
REAL ( dbl ), DIMENSION ( : ), INTENT ( OUT ) :: y
REAL ( dbl ), DIMENSION ( :, : ), INTENT ( OUT ) :: p
! Locals
INTEGER :: i, j, iochannel = 99
!------------------------------------------------------------------------------
! read from the restart file
OPEN ( UNIT = iochannel, FILE = 'amoeba.rst', STATUS = 'UNKNOWN' )
DO i = 1, SIZE ( y )
READ ( iochannel, FMT = * ) y ( i )
END DO
DO i = 1, SIZE ( y )
READ ( iochannel, FMT = * ) ( p ( i, j ), j = 1, SIZE ( p, 2 ) )
END DO
CLOSE ( UNIT = iochannel )
END SUBROUTINE fit_read_restart
!******************************************************************************
!
END MODULE fit_utilities
!------------------------------------------------------------------------------!

133
src/force_control.F
View file

@ -26,6 +26,7 @@ MODULE force_control
USE dft_structure_control, ONLY : dft_force_control
USE dft_types, ONLY : dft_control_type
USE fist_force, ONLY : fist_force_control => force_control
USE pol_force, ONLY : pol_force_control => force_control
USE global_types, ONLY : global_environment_type
USE tbmd_force, ONLY : tbmd_force_control => force_control
USE kinds, ONLY : dbl
@ -41,7 +42,6 @@ MODULE force_control
INTERFACE force
MODULE PROCEDURE force_1, force_n
END INTERFACE
!!*****
!******************************************************************************
@ -66,13 +66,13 @@ CONTAINS
!! SOURCE
!******************************************************************************
SUBROUTINE force_n ( struc, inter, thermo, simpar, ewald_param, box_change, &
SUBROUTINE force_1 ( struc, inter, thermo, simpar, ewald_param, box_change, &
globenv, dft_control, kp )
IMPLICIT NONE
! Arguments
TYPE ( structure_type ),DIMENSION ( : ), INTENT ( INOUT ) :: struc
TYPE ( structure_type ), INTENT ( INOUT ) :: struc
TYPE ( thermodynamic_type ), INTENT ( INOUT ) :: thermo
TYPE ( simulation_parameters_type ), INTENT ( IN ) :: simpar
TYPE ( ewald_parameters_type ), INTENT ( INOUT ) :: ewald_param
@ -81,7 +81,62 @@ SUBROUTINE force_n ( struc, inter, thermo, simpar, ewald_param, box_change, &
TYPE ( global_environment_type), INTENT ( IN ) :: globenv
TYPE ( dft_control_type ), INTENT ( IN ) :: dft_control
TYPE ( kpoint_type ), INTENT ( IN ) :: kp
!------------------------------------------------------------------------------
thermo % ptens_aval = .FALSE.
SELECT CASE ( simpar % program )
CASE DEFAULT
CALL stop_program ( 'force', 'not implemented' )
CASE ( 'FIST' )
CALL fist_force_control ( struc % molecule, struc % pnode, struc % part, &
struc % box, struc % box_ref, thermo, inter%potparm, ewald_param, &
box_change, globenv, struc%ll_data (1) )
CASE ( 'POL' )
CALL pol_force_control ( struc % molecule, struc % pnode, struc % part, &
struc % box, struc % box_ref, struc % drho_basis_info, &
struc % rho0_basis_info, struc % coef_pos (1), struc % coef_vel (1), &
struc % coef_force (1), thermo, inter%potparm, ewald_param, &
simpar % ensemble, globenv, struc%ll_data )
CASE ( 'TBMD' )
CALL tbmd_force_control ( struc % molecule, struc % pnode, struc % part, &
struc % box, struc % box_ref, thermo, inter%potparm, ewald_param, &
box_change, globenv , struc%ll_data (1) )
CASE ( 'WAVE' )
CALL dft_force_control ( struc % molecule, struc % pnode, struc % part, &
struc % box, thermo, inter % potparm, ewald_param, &
box_change, struc % ll_data (1), dft_control, kp, &
dftpar = globenv, wavepar = globenv )
END SELECT
END SUBROUTINE force_1
!!*****
!******************************************************************************
SUBROUTINE force_n ( struc, inter, thermo, simpar, ewald_param, box_change, &
globenv, dft_control, kp )
IMPLICIT NONE
! Arguments
TYPE ( structure_type ), DIMENSION ( : ), INTENT ( INOUT ) :: struc
TYPE ( thermodynamic_type ), INTENT ( INOUT ) :: thermo
TYPE ( simulation_parameters_type ), INTENT ( IN ) :: simpar
TYPE ( ewald_parameters_type ), INTENT ( INOUT ) :: ewald_param
LOGICAL, INTENT ( IN ) :: box_change
TYPE ( interaction_type ), INTENT ( IN ) :: inter
TYPE ( global_environment_type), INTENT ( IN ) :: globenv
TYPE ( dft_control_type ), INTENT ( IN ) :: dft_control
TYPE ( kpoint_type ), INTENT ( IN ) :: kp
!------------------------------------------------------------------------------
! Local
INTEGER :: i, n
REAL (dbl) :: e_pot
@ -95,26 +150,24 @@ SUBROUTINE force_n ( struc, inter, thermo, simpar, ewald_param, box_change, &
e_pot = 0.0_dbl
DO i = 1, n
SELECT CASE ( simpar % program )
CASE DEFAULT
CALL stop_program ( 'force', 'not implemented' )
CASE ( 'FIST' )
CALL fist_force_control ( struc(i) % molecule, struc(i) % pnode, &
struc(i) % part, struc(i) % box, struc(i) % box_ref, thermo, &
inter%potparm, ewald_param, box_change, globenv, struc(i) % ll_data )
CALL fist_force_control ( struc ( i ) % molecule, struc ( i ) % pnode, &
struc ( i ) % part, struc ( i ) % box, struc ( i ) % box_ref, &
thermo, inter%potparm, ewald_param, box_change, globenv, &
struc ( i ) % ll_data ( 1 ) )
CASE ( 'POL' )
CALL stop_program ( 'force', 'POL PIMD not implemented' )
CASE ( 'TBMD' )
CALL tbmd_force_control ( struc(i) % molecule, struc(i) % pnode, &
struc(i) % part, struc(i) % box, struc(i) % box_ref, thermo, &
inter%potparm, ewald_param, box_change, globenv, struc(i) % ll_data )
CALL stop_program ( 'force', 'TBMD PIMD not implemented' )
CASE ( 'WAVE' )
CALL dft_force_control ( struc(i) % molecule, struc(i) % pnode, &
struc(i) % part, struc(i) % box, thermo, inter % potparm, &
ewald_param, box_change, struc(i) % ll_data, dft_control, kp, &
dftpar = globenv, wavepar = globenv )
CALL stop_program ( 'force', 'WAVE PIMD not implemented' )
END SELECT
@ -125,56 +178,6 @@ SUBROUTINE force_n ( struc, inter, thermo, simpar, ewald_param, box_change, &
thermo % pot = e_pot
END SUBROUTINE force_n
!******************************************************************************
SUBROUTINE force_1 ( struc, inter, thermo, simpar, ewald_param, box_change, &
globenv, dft_control, kp )
IMPLICIT NONE
! Arguments
TYPE ( structure_type ), INTENT ( INOUT ) :: struc
TYPE ( thermodynamic_type ), INTENT ( INOUT ) :: thermo
TYPE ( simulation_parameters_type ), INTENT ( IN ) :: simpar
TYPE ( ewald_parameters_type ), INTENT ( INOUT ) :: ewald_param
LOGICAL, INTENT ( IN ) :: box_change
TYPE ( interaction_type ), INTENT ( IN ) :: inter
TYPE ( global_environment_type), INTENT ( IN ) :: globenv
TYPE ( dft_control_type ), INTENT ( IN ) :: dft_control
TYPE ( kpoint_type ), INTENT ( IN ) :: kp
!------------------------------------------------------------------------------
thermo % ptens_aval = .FALSE.
SELECT CASE ( simpar % program )
CASE DEFAULT
CALL stop_program ( 'force', 'not implemented' )
CASE ( 'FIST' )
CALL fist_force_control ( struc % molecule, struc % pnode, &
struc % part, struc % box, struc % box_ref, thermo, &
inter%potparm, ewald_param, box_change, globenv, struc % ll_data )
CASE ( 'TBMD' )
CALL tbmd_force_control ( struc % molecule, struc % pnode, &
struc % part, struc % box, struc % box_ref, thermo, &
inter%potparm, ewald_param, box_change, globenv, struc % ll_data )
CASE ( 'WAVE' )
CALL dft_force_control ( struc % molecule, struc % pnode, &
struc % part, struc % box, thermo, inter % potparm, &
ewald_param, box_change, struc % ll_data, dft_control, kp, &
dftpar = globenv, wavepar = globenv )
END SELECT
END SUBROUTINE force_1
!!*****
!******************************************************************************
END MODULE force_control
!******************************************************************************

46
src/function_amoeba.F Normal file
View file

@ -0,0 +1,46 @@
!------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 1999 MPI fuer Festkoerperforschung, Stuttgart !
!------------------------------------------------------------------------------!
!*****************************************************************
!!
!! NAME
!! function_amoeba
!!
!*****************************************************************
!------------------------------------------------------------------------------!
FUNCTION func(x,struc,potparm,thermo,simpar,fit_param,box_ref, &
finfo,param,pstat,ewald_param, globenv)
USE kinds, ONLY : dbl
USE structure_types, ONLY : structure_type
USE pair_potential, ONLY : potentialparm_type
USE particle_types, ONLY : particle_prop_type
USE md, ONLY : simulation_parameters_type, thermodynamic_type
USE simulation_cell, ONLY : cell_type
USE fit_parameters_types, ONLY : fit_parameter_type, linklist_parameters_type, &
fit_info_type
USE ewald_parameters_types, ONLY : ewald_parameters_type
USE global_types, ONLY : global_environment_type
USE fit_module, ONLY : local
IMPLICIT NONE
TYPE (structure_type), INTENT (INOUT) :: struc
TYPE (potentialparm_type), DIMENSION (:,:), INTENT (INOUT) :: potparm
TYPE (thermodynamic_type), INTENT (INOUT) :: thermo
TYPE (simulation_parameters_type), INTENT (IN) :: simpar
REAL (dbl), INTENT (IN), DIMENSION (:) :: x
TYPE (fit_parameter_type), INTENT(INOUT) :: fit_param
TYPE (cell_type) , POINTER, DIMENSION (:) :: box_ref
TYPE (fit_info_type), DIMENSION (:,:), POINTER :: finfo
TYPE (linklist_parameters_type ), INTENT(INOUT), DIMENSION ( :, : ) :: param
TYPE (particle_prop_type ), DIMENSION ( : ), POINTER :: pstat
TYPE ( ewald_parameters_type ), INTENT ( INOUT ) :: ewald_param
TYPE ( global_environment_type), INTENT ( IN ) :: globenv
REAL (dbl) :: func
func = local(x,struc,potparm,thermo,simpar,fit_param,box_ref, &
finfo,param,pstat,ewald_param, globenv)
END FUNCTION func

10
src/hartree.F
View file

@ -56,8 +56,8 @@ CONTAINS
SUBROUTINE calculate_hartree ( density, ehartree, vhartree, dvhartree, h_stress )
! Arguments
TYPE ( coeff_type ), INTENT ( INOUT ) :: density
REAL ( dbl ), INTENT ( OUT ) :: ehartree
TYPE ( coeff_type ), INTENT ( IN ) :: density
REAL ( dbl ), INTENT ( OUT ), OPTIONAL :: ehartree
TYPE ( coeff_type ), TARGET, INTENT ( OUT ), OPTIONAL :: vhartree
TYPE ( coeff_type ), DIMENSION ( 3 ), INTENT ( OUT ), OPTIONAL :: dvhartree
REAL ( dbl ), DIMENSION ( 3, 3 ), INTENT ( OUT ), OPTIONAL :: h_stress
@ -114,9 +114,9 @@ SUBROUTINE poisson_solver ( density, ehartree, vhartree, dvhartree, h_stress )
IMPLICIT NONE
! Arguments
TYPE ( coeff_type ), INTENT ( INOUT ) :: density
REAL ( dbl ), INTENT ( OUT ) :: ehartree
TYPE ( coeff_type ), INTENT ( OUT ) :: vhartree
TYPE ( coeff_type ), INTENT ( IN ) :: density
REAL ( dbl ), INTENT ( OUT ), OPTIONAL :: ehartree
TYPE ( coeff_type ), INTENT ( OUT ), OPTIONAL :: vhartree
TYPE ( coeff_type ), DIMENSION ( 3 ), INTENT ( OUT ), OPTIONAL :: dvhartree
REAL ( dbl ), DIMENSION ( 3, 3 ), INTENT ( OUT ), OPTIONAL :: h_stress

215
src/initialize_extended_types.F
View file

@ -14,27 +14,30 @@
!!
!! MODIFICATION HISTORY
!! cjm, FEB 20 2001: added subroutine assign_extended_parameters
!!
!! cjm, MAY 03 2001: reorganized and added separtate routines for
!! nhc_part, nhc_baro, nhc_ao, npt
!! SOURCE
!******************************************************************************
MODULE initialize_extended_types
USE atoms_input, ONLY : system_type
USE coefficient_types, ONLY : coeff_type
USE global_types, ONLY : global_environment_type
USE kinds, ONLY : dbl
USE md, ONLY : simulation_parameters_type
USE molecule_types, ONLY : molecule_structure_type, molecule_type
USE message_passing, ONLY : mp_sum
USE nhc, ONLY : nhc_control_particles, nhc_control_barostat
USE nose, ONLY : initial_nhc, lnhc_parameters_type, initial_npt, &
USE nhc, ONLY : nhc_particle_setup, nhc_barostat_setup, nhc_pol_setup
USE nose, ONLY : initialize_nhc, lnhc_parameters_type, initialize_npt, &
npt_info_type
USE simulation_cell, ONLY : cell_type
USE termination, ONLY : stop_program
IMPLICIT NONE
PUBLIC :: initialize_extended_system, assign_extended_parameters
PUBLIC :: initialize_nhc_part, initialize_nhc_baro, initialize_npt_type, &
initialize_nhc_pol, assign_extended_parameters
PRIVATE
!!*****
@ -43,10 +46,10 @@ MODULE initialize_extended_types
CONTAINS
!******************************************************************************
!!****** initialize_extended_types/initialize_extended_system [1.0] *
!!****** initialize_extended_types/initialize_npt_type [1.0] *
!!
!! NAME
!! initialize_extended_types
!! initialize_npt_type
!!
!! FUNCTION
!!
@ -54,61 +57,25 @@ CONTAINS
!! CJM
!!
!! MODIFICATION HISTORY
!! cjm FEB-20-2001: Only initializes when restart type
!! is init, pos, posvel
!!
!!
!! SOURCE
!******************************************************************************
SUBROUTINE initialize_extended_system ( box, simpar, molecule, mol_setup, &
globenv, nhc_part, nhc_baro, npt_info )
SUBROUTINE initialize_npt_type ( simpar, globenv, npt_info, box )
IMPLICIT NONE
! Arguments
TYPE ( simulation_parameters_type ), INTENT ( INOUT ) :: simpar
TYPE ( molecule_structure_type ), DIMENSION ( : ), INTENT ( INOUT ) :: &
molecule
TYPE ( molecule_type ), DIMENSION ( : ), INTENT ( INOUT ) :: mol_setup
TYPE ( cell_type ), INTENT ( IN ) :: box
TYPE ( global_environment_type ), INTENT ( INOUT ) :: globenv
TYPE ( lnhc_parameters_type ), INTENT ( INOUT ) :: nhc_part
TYPE ( lnhc_parameters_type ), INTENT ( INOUT ) :: nhc_baro
TYPE ( npt_info_type ), DIMENSION ( :, : ), POINTER :: npt_info
TYPE ( cell_type ), INTENT (inout) :: box
! locals
INTEGER :: isos
!------------------------------------------------------------------------------
! fire up the thermostats, if not NVE
CALL nhc_control_particles ( simpar, molecule, mol_setup, nhc_part, globenv )
! fire up the thermostats, if NPT
CALL nhc_control_barostat ( simpar, nhc_baro, globenv )
IF ( simpar % read_type /= 'ALL' ) THEN
! Initializing thermostat forces and velocities for the Nose-Hoover
! Chain variables
CALL initial_nhc ( nhc_part, simpar % tau_nhc, &
simpar % temp_ext, globenv )
CALL initial_nhc ( nhc_baro, simpar % tau_nhc, &
simpar % temp_ext, globenv )
END IF
! compute the total number of degrees of freedom for temperature
! Initialize simpar%nfree to dim (because at this point no
! assumptions have been made about thermostatting options
! and COM is ALWAYS zeroed in this code. If there are
! multiple thermostats etc...simpar%nfree will be set
! to zero in nhc_control.f)
simpar % nfree = 6 - SUM ( box % perd ( 1:3 ) )
! first allocating the npt_info_type
SELECT CASE ( simpar % ensemble )
CASE ( 'NPT_I')
@ -127,13 +94,161 @@ SUBROUTINE initialize_extended_system ( box, simpar, molecule, mol_setup, &
'failed to allocate npt_info')
END SELECT
CALL get_nfree ( simpar, molecule, globenv )
IF ( simpar % read_type /= 'ALL' ) &
CALL initial_npt ( npt_info, simpar % tau_cell, simpar % temp_ext, &
CALL initialize_npt ( npt_info, simpar % tau_cell, simpar % temp_ext, &
box % deth, simpar % nfree, simpar % ensemble, globenv )
END SUBROUTINE initialize_extended_system
END SUBROUTINE initialize_npt_type
!******************************************************************************
!!****** initialize_extended_types/initialize_nhc_pol [1.0] *
!!
!! NAME
!! initialize_nhc_pol
!!
!! FUNCTION
!!
!! AUTHOR
!! CJM
!!
!! MODIFICATION HISTORY
!!
!!
!! SOURCE
!******************************************************************************
SUBROUTINE initialize_nhc_pol ( simpar, coeff, globenv, nhc )
IMPLICIT NONE
! Arguments
TYPE ( coeff_type ), INTENT ( IN ) :: coeff
TYPE ( simulation_parameters_type ), INTENT ( INOUT ) :: simpar
TYPE ( global_environment_type ), INTENT ( INOUT ) :: globenv
TYPE ( lnhc_parameters_type ), INTENT ( INOUT ) :: nhc
! locals
INTEGER :: isos
!------------------------------------------------------------------------------
! fire up the thermostats
CALL nhc_pol_setup ( simpar, coeff, nhc, globenv )
IF ( simpar % pol_read_type /= 'ALL' ) THEN
! Initializing thermostat forces and velocities for the Nose-Hoover
! Chain variables
IF ( nhc % nhc_len /= 0 ) CALL initialize_nhc ( nhc, simpar % pol_tau_nhc, &
simpar % pol_temp, globenv )
END IF
END SUBROUTINE initialize_nhc_pol
!******************************************************************************
!!****** initialize_extended_types/initialize_nhc_baro [1.0] *
!!
!! NAME
!! initialize_nhc_baro
!!
!! FUNCTION
!!
!! AUTHOR
!! CJM
!!
!! MODIFICATION HISTORY
!!
!!
!! SOURCE
!******************************************************************************
SUBROUTINE initialize_nhc_baro ( simpar, globenv, nhc )
IMPLICIT NONE
! Arguments
TYPE ( simulation_parameters_type ), INTENT ( INOUT ) :: simpar
TYPE ( global_environment_type ), INTENT ( INOUT ) :: globenv
TYPE ( lnhc_parameters_type ), INTENT ( INOUT ) :: nhc
! locals
INTEGER :: isos
!------------------------------------------------------------------------------
! fire up the thermostats, if NPT
CALL nhc_barostat_setup ( simpar, nhc, globenv )
IF ( simpar % read_type /= 'ALL' ) THEN
! Initializing thermostat forces and velocities for the Nose-Hoover
! Chain variables
IF ( nhc % nhc_len /= 0 ) &
CALL initialize_nhc ( nhc, simpar % tau_nhc, &
simpar % temp_ext, globenv )
END IF
END SUBROUTINE initialize_nhc_baro
!******************************************************************************
!!****** initialize_extended_types/initialize_nhc_part [1.0] *
!!
!! NAME
!! initialize_nhc_part
!!
!! FUNCTION
!!
!! AUTHOR
!! CJM
!!
!! MODIFICATION HISTORY
!!
!! SOURCE
!******************************************************************************
SUBROUTINE initialize_nhc_part ( box, simpar, molecule, mol_setup, &
globenv, nhc )
IMPLICIT NONE
! Arguments
TYPE ( simulation_parameters_type ), INTENT ( INOUT ) :: simpar
TYPE ( molecule_structure_type ), DIMENSION ( : ), INTENT ( INOUT ) :: &
molecule
TYPE ( molecule_type ), DIMENSION ( : ), INTENT ( INOUT ) :: mol_setup
TYPE ( cell_type ), INTENT ( IN ) :: box
TYPE ( global_environment_type ), INTENT ( INOUT ) :: globenv
TYPE ( lnhc_parameters_type ), INTENT ( INOUT ) :: nhc
! locals
INTEGER :: isos
!------------------------------------------------------------------------------
! Initialize simpar%nfree to dim (because at this point no
! assumptions have been made about thermostatting options
! and COM is ALWAYS zeroed in this code. If there are
! multiple thermostats etc...simpar % nfree will be set
! to zero in nhc_particle_setup.f)
simpar % nfree = 6 - SUM ( box % perd ( 1:3 ) )
! fire up the thermostats, if not NVE
CALL nhc_particle_setup ( simpar, molecule, mol_setup, nhc, globenv )
IF ( simpar % read_type /= 'ALL' ) THEN
! Initializing thermostat forces and velocities for the Nose-Hoover
! Chain variables
IF ( nhc % nhc_len /= 0 ) &
CALL initialize_nhc ( nhc, simpar % tau_nhc, &
simpar % temp_ext, globenv )
ENDIF
! compute the total number of degrees of freedom for temperature
CALL get_nfree ( simpar, molecule, globenv )
END SUBROUTINE initialize_nhc_part
!!*****
!******************************************************************************
@ -276,8 +391,6 @@ SUBROUTINE assign_extended_parameters ( deth, ensemble, ainp, npt, &
IF ( ensemble == 'NPT_I' ) npt ( :, : ) % eps = LOG ( deth ) / 3.0_dbl
END SUBROUTINE assign_extended_parameters
!!*****

4
src/initialize_molecule_types.F
View file

@ -101,6 +101,7 @@ SUBROUTINE initialize_molecule_type ( mol_setup, intra_param, &
END DO
DO iat = 1, size(pnode)
pnode(iat) %p => part(iat)
part(iat)%node = iat
END DO
ELSE
@ -138,7 +139,8 @@ SUBROUTINE initialize_molecule_type ( mol_setup, intra_param, &
atombase = molecule(imol) % atombase - 1
DO j = 1, molecule(imol) % natoms_mol
k = k + 1
pnode(k) % p => part(atombase+j)
pnode ( k ) % p => part ( atombase + j )
part ( atombase + j ) % node = k
END DO
END DO
END IF

10
src/initialize_particle_types.F
View file

@ -78,11 +78,15 @@ SUBROUTINE initialize_particle_type ( atom_names, simpar, mol_setup, &
part(iat) % r(3) = ainp % c(3,iat)
END DO
part ( : ) % v ( 1 ) = 0._dbl
part ( : ) % v ( 2 ) = 0._dbl
part ( : ) % v ( 3 ) = 0._dbl
IF ( ( simpar%read_type=='POSVEL' ) .OR. ( simpar%read_type=='ALL' ) ) THEN
DO iat = 1, natoms
part(iat) %v(1) = ainp%v(1,iat)
part(iat) %v(2) = ainp%v(2,iat)
part(iat) %v(3) = ainp%v(3,iat)
part ( iat ) % v ( 1 ) = ainp % v ( 1, iat )
part ( iat ) % v ( 2 ) = ainp % v ( 2, iat )
part ( iat ) % v ( 3 ) = ainp % v ( 3, iat )
END DO
END IF

1288
src/integrator.F
View file

File diff suppressed because it is too large Load diff

39
src/lib/fast.F
View file

@ -24,19 +24,46 @@ SUBROUTINE rankup ( n, za, cmat, zb, ex, ey, ez, scr )
n3 = n2 * n ( 3 )
scr ( 1:n2 ) = CMPLX ( 0._dbl, KIND = dbl )
#if defined (__sp_lib)
CALL cgeru ( n ( 1 ), n ( 2 ), zb, ex, 1, ey, 1, scr, n ( 1 ) )
CALL cscal ( n3, za, cmat, 1 )
CALL cgeru ( n2, n ( 3 ), cone, scr, 1, ez, 1, cmat, n2 )
CALL CGERU ( n ( 1 ), n ( 2 ), zb, ex, 1, ey, 1, scr, n ( 1 ) )
CALL CSCAL ( n3, za, cmat, 1 )
CALL CGERU ( n2, n ( 3 ), cone, scr, 1, ez, 1, cmat, n2 )
#else
CALL zgeru ( n ( 1 ), n ( 2 ), zb, ex, 1, ey, 1, scr, n ( 1 ) )
CALL zscal ( n3, za, cmat, 1 )
CALL zgeru ( n2, n ( 3 ), cone, scr, 1, ez, 1, cmat, n2 )
CALL ZGERU ( n ( 1 ), n ( 2 ), zb, ex, 1, ey, 1, scr, n ( 1 ) )
CALL ZSCAL ( n3, za, cmat, 1 )
CALL ZGERU ( n2, n ( 3 ), cone, scr, 1, ez, 1, cmat, n2 )
#endif
END SUBROUTINE rankup
!-----------------------------------------------------------------------------!
SUBROUTINE vc_x_vc ( n, cvec2, cvec )
!
! cden(i) <- cden(i) * cden2(i)
!
IMPLICIT NONE
INTEGER, PARAMETER :: dbl = SELECTED_REAL_KIND ( 14, 200 )
INTEGER, INTENT ( IN ) :: n
REAL ( dbl ), INTENT ( IN ) :: cvec2 (2,*)
REAL ( dbl ), INTENT ( INOUT ) :: cvec (2,*)
INTEGER :: i
REAL ( dbl ) :: re_c, im_c, re_c2, im_c2
DO i = 1, n
re_c = cvec ( 1, i )
im_c = cvec ( 2, i )
re_c2 = cvec2 ( 1, i )
im_c2 = cvec2 ( 2, i )
cvec ( 1, i ) = re_c * re_c2 - im_c * im_c2
cvec ( 2, i ) = im_c * re_c2 + re_c * im_c2
END DO
END SUBROUTINE vc_x_vc
!-----------------------------------------------------------------------------!
SUBROUTINE vr_x_vc ( n, rvec, cvec )
!
! cden(i) <- cden(i) * rden(i)

4
src/lib/mltfftsg.F
View file

@ -49,10 +49,10 @@ SUBROUTINE mltfftsg ( transa, transb, a, ldax, lday, b, ldbx, ldby, &
ENDIF
IF ( TSCAL ) THEN
IF ( LOT == NFFT ) THEN
CALL dscal ( 2 * LOT * N, SCALE, Z ( 1, 1 ), 1 )
CALL DSCAL ( 2 * LOT * N, SCALE, Z ( 1, 1 ), 1 )
ELSE
DO I = 1, N
CALL dscal ( 2 * NFFT, SCALE, Z ( LOT * ( I - 1 ) + 1, 1 ), 1 )
CALL DSCAL ( 2 * NFFT, SCALE, Z ( LOT * ( I - 1 ) + 1, 1 ), 1 )
END DO
END IF
END IF

10
src/library_tests.F
View file

@ -381,7 +381,7 @@ SUBROUTINE matmul_test ( globenv )
tstart = m_cputime ( )
DO j = 1, ntim
CALL dgemm ( "N", "N", len, len, len, 1._dbl, ma, len, mb, len, 0._dbl, mc, len )
CALL DGEMM ( "N", "N", len, len, len, 1._dbl, ma, len, mb, len, 0._dbl, mc, len )
END DO
tend = m_cputime ( )
t = tend - tstart
@ -394,8 +394,8 @@ SUBROUTINE matmul_test ( globenv )
tstart = m_cputime ( )
DO j = 1, ntim
CALL dgemm ( "N", "N", len, len, len, 1._dbl, ma, len, mb, len, 0._dbl, mc, len )
CALL dcopy ( len * len , mc, 1, ma, 1)
CALL DGEMM ( "N", "N", len, len, len, 1._dbl, ma, len, mb, len, 0._dbl, mc, len )
CALL DCOPY ( len * len , mc, 1, ma, 1)
END DO
tend = m_cputime ( )
t = tend - tstart
@ -408,7 +408,7 @@ SUBROUTINE matmul_test ( globenv )
tstart = m_cputime ( )
DO j = 1, ntim
CALL dgemm ( "N", "T", len, len, len, 1._dbl, ma, len, mb, len, 0._dbl, mc, len )
CALL DGEMM ( "N", "T", len, len, len, 1._dbl, ma, len, mb, len, 0._dbl, mc, len )
END DO
tend = m_cputime ( )
t = tend - tstart
@ -421,7 +421,7 @@ SUBROUTINE matmul_test ( globenv )
tstart = m_cputime ( )
DO j = 1, ntim
CALL dgemm ( "T", "N", len, len, len, 1._dbl, ma, len, mb, len, 0._dbl, mc, len )
CALL DGEMM ( "T", "N", len, len, len, 1._dbl, ma, len, mb, len, 0._dbl, mc, len )
END DO
tend = m_cputime ( )
t = tend - tstart

12
src/linklist_cell_list.F
View file

@ -46,7 +46,7 @@ CONTAINS
! used in the construction of all neighbor lists
SUBROUTINE cell_list ( natom_types, n_cell, rcut, tag, pnode, part, &
box, cell_ll, start_list, neighbor, globenv )
box, cell_ll, start_list, neighbor, globenv, list_type )
IMPLICIT NONE
@ -62,6 +62,7 @@ SUBROUTINE cell_list ( natom_types, n_cell, rcut, tag, pnode, part, &
INTEGER, INTENT ( IN ), DIMENSION ( : ) :: n_cell
INTEGER, INTENT ( IN ) :: natom_types
CHARACTER ( LEN = * ), INTENT ( IN ) :: tag
INTEGER, INTENT(IN) :: list_type
! Locals
TYPE ( cell_atoms_linklist_type ), POINTER :: loc_ll
@ -216,6 +217,8 @@ SUBROUTINE cell_list ( natom_types, n_cell, rcut, tag, pnode, part, &
DO i = 1, natoms
itype = part(i) % prop % ptype
IF (itype > natom_types) CYCLE
! pbc so that particles outside the central box are also mapped
! onto proper cells
r(1) = part(i) % r(1)
@ -258,7 +261,12 @@ SUBROUTINE cell_list ( natom_types, n_cell, rcut, tag, pnode, part, &
END IF
END DO
loc_ll % atom = i
part(i) % cell => cell_ll(bin(1),bin(2),bin(3))
SELECT CASE (list_type)
CASE (1)
part(i) % cell => cell_ll(bin(1),bin(2),bin(3))
CASE (2)
part(i) % cell_ol => cell_ll(bin(1),bin(2),bin(3))
END SELECT
END DO
! linked list that forms a list of all neighboring cells (done only once)

32
src/linklist_control.F
View file

@ -21,6 +21,7 @@
MODULE linklist_control
USE atomic_kinds, ONLY : kind_info_type
USE global_types, ONLY : global_environment_type
USE kinds, ONLY : dbl
USE linklist_cell_list, ONLY : cell_list
@ -29,6 +30,7 @@ MODULE linklist_control
neighbor_cells_type
USE linklist_utilities, ONLY : startlist_images_type
USE linklist_verlet_list, ONLY : verlet_list
USE linklist_overlap_list, ONLY : overlap_list
USE molecule_types, ONLY : particle_node_type,linklist_neighbor, &
linklist_exclusion, linklist_images
USE message_passing, ONLY : mp_max
@ -124,6 +126,8 @@ SUBROUTINE set_ll_parm ( internal_data, globenv_in, vskin, natype, rcut, &
internal_data%list_type = 0
CASE ( 'NONBOND')
internal_data%list_type = 1
CASE ( 'OVERLAP')
internal_data%list_type = 2
END SELECT
END IF
@ -151,7 +155,7 @@ END SUBROUTINE get_ll_parm
! if any particle has moved by 0.5*skin_length from the previous
! list_update, then the list routine is called.
!
SUBROUTINE list_control ( internal_data, pnode, part, box )
SUBROUTINE list_control ( internal_data, pnode, part, box, ki )
IMPLICIT NONE
@ -160,6 +164,7 @@ SUBROUTINE list_control ( internal_data, pnode, part, box )
TYPE ( particle_type ), DIMENSION ( : ), INTENT ( INOUT ) :: part
TYPE ( particle_node_type ), DIMENSION ( : ), INTENT ( INOUT ) :: pnode
TYPE ( cell_type ), INTENT ( IN ) :: box
TYPE (kind_info_type), INTENT (IN), DIMENSION (:), OPTIONAL :: ki
! Locals
REAL ( dbl ), DIMENSION (3) :: disp_vec
@ -215,14 +220,35 @@ SUBROUTINE list_control ( internal_data, pnode, part, box )
CALL cell_list(internal_data%natom_types,internal_data%pp_cell,&
internal_data%rlist_cut,string,pnode,part, &
box,internal_data%pp_cell_ll,internal_data%pp_startlist_cell,&
internal_data%pp_neighbor, internal_data%globenv_ll)
internal_data%pp_neighbor, internal_data%globenv_ll, &
internal_data%list_type)
CALL verlet_list(internal_data%natom_types,internal_data%pp_images,&
internal_data%pp_cell,pnode,part,box, &
internal_data%verlet_skin,internal_data%rlist_cut,&
internal_data%rlist_cutsq,internal_data%pp_neighbor,&
internal_data%pp_cell_ll,internal_data%pp_startlist_nl,&
internal_data%pp_startlist_im, internal_data%globenv_ll)
internal_data%pp_startlist_im, internal_data%globenv_ll, &
internal_data%list_type)
CASE ( 2 )
string = ' cell_index: overlap neighbor list'
CALL cell_list(internal_data%natom_types,internal_data%pp_cell,&
internal_data%rlist_cut,string,pnode,part, &
box,internal_data%pp_cell_ll,internal_data%pp_startlist_cell,&
internal_data%pp_neighbor, internal_data%globenv_ll, &
internal_data%list_type)
CALL overlap_list(internal_data%pp_images,&
internal_data%pp_cell,pnode,part,box, ki, &
internal_data%verlet_skin,internal_data%rlist_cut,&
internal_data%rlist_cutsq,internal_data%pp_neighbor,&
internal_data%pp_cell_ll,internal_data%pp_startlist_nl,&
internal_data%pp_startlist_im, internal_data%globenv_ll, &
internal_data%list_type)
END SELECT
iw = internal_data%globenv_ll % scr

370
src/linklist_overlap_list.F Normal file
View file

@ -0,0 +1,370 @@
!-----------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/linklist_overlap_list [1.0] *
!!
!! NAME
!! linklist_overlap_list
!!
!! FUNCTION
!!
!! AUTHOR
!! CJM
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE linklist_overlap_list
USE kinds, ONLY : dbl
USE atomic_kinds, ONLY : kind_info_type
USE global_types, ONLY : global_environment_type
USE linklist_cell_types, ONLY : cell_neighborlist_type, &
neighbor_cells_type, cell_atoms_linklist_type
USE linklist_utilities, ONLY : get_match, update_verlet_list, &
number_of_image_cells, startlist_images_type
USE molecule_types, ONLY : particle_node_type, linklist_neighbor, &
linklist_exclusion
USE message_passing, ONLY : mp_max, mp_min, mp_sum
USE particle_types, ONLY : particle_type
USE simulation_cell, ONLY : cell_type
USE termination, ONLY : stop_memory
USE timings, ONLY : timeset, timestop
USE util, ONLY : include_list, sort
IMPLICIT NONE
PRIVATE
PUBLIC :: overlap_list
!!*****
!******************************************************************************
CONTAINS
!******************************************************************************
! to get the Verlet neighbor list for nonbonded interactions
SUBROUTINE overlap_list(n_images,n_cell,pnode,part,box, &
ki,verlet_skin,rlist_cut,rlist_cutsq,neighbor, &
cell_ll,start_list_nl,start_list_im, globenv,list_type)
IMPLICIT NONE
!
! neighbor_list lists all the neighbors of all atoms
!
! Arguments
TYPE ( particle_node_type ), DIMENSION ( : ), INTENT ( INOUT ), TARGET :: &
pnode
TYPE ( particle_type ), DIMENSION ( : ), INTENT ( IN ), TARGET :: part
TYPE ( cell_type ), INTENT ( IN ) :: box
TYPE ( kind_info_type), INTENT (IN), DIMENSION (:) :: ki
TYPE ( linklist_neighbor), POINTER :: start_list_nl
TYPE ( startlist_images_type ), POINTER :: start_list_im ( : )
TYPE ( cell_neighborlist_type ), DIMENSION ( :, :, : ), POINTER :: cell_ll
TYPE ( neighbor_cells_type ), DIMENSION ( :, : ), POINTER :: neighbor
TYPE ( global_environment_type ), INTENT ( IN ) :: globenv
REAL ( dbl ), DIMENSION ( :, : ), INTENT ( IN ) :: rlist_cut, rlist_cutsq
REAL ( dbl ), INTENT ( IN ) :: verlet_skin
INTEGER, DIMENSION ( :, :, : ), POINTER :: n_images
INTEGER, DIMENSION ( : ), INTENT ( IN ) :: n_cell
INTEGER, INTENT (IN) :: list_type
! Locals
TYPE (cell_atoms_linklist_type ), POINTER :: ll
INTEGER :: i, ii, j, natoms, nnodes, handle, ipair
INTEGER :: iatom, jatom, iexclude, isos
INTEGER :: nkind, ikind, jkind, nat
INTEGER :: nloc, nsum, nmax, nmin, nexcl, ncells
INTEGER :: k, ipart, jpart, npairs
INTEGER, DIMENSION (3) :: index
INTEGER, PARAMETER :: els = 50
INTEGER, DIMENSION ( : ), POINTER :: elist, ind
REAL ( dbl ) :: rijsq, fis
REAL ( dbl ), DIMENSION (3) :: perd
REAL ( dbl ), DIMENSION (3,3) :: hmat, h_inv
LOGICAL :: first_time
!
! the following lines are used for the link_list of
! neighbors
!
TYPE (linklist_neighbor), POINTER :: current_neighbor
TYPE (linklist_exclusion), POINTER :: ll_exclude
!------------------------------------------------------------------------------
! initialize values
natoms = 0
nat = 0
!
! if we are entering this subroutine for the first time
! create the start of the link_list of neighbors.
! start_list and first_time are both of type save
! so this if statement will not be run or needed
! after the first call to the subroutine
!
CALL timeset ( 'VERLET_LIST', 'I', 'KIs', handle )
fis = 0.0_dbl
ALLOCATE (elist(els),STAT=isos)
IF ( isos /= 0 ) CALL stop_memory ( 'overlap_list', 'elist', els )
ALLOCATE (ind(els),STAT=isos)
IF ( isos /= 0 ) CALL stop_memory ( 'overlap_list', 'ind', els )
! local copies of box parameters
hmat = box%hmat
h_inv = box%h_inv
perd = box%perd
nnodes = size(pnode)
nkind = size(KI)
! count number of atoms with basis functions
DO i = 1, nkind
natoms = natoms + size ( ki(i)%atom_list )
END DO
first_time = .NOT. ASSOCIATED ( n_images )
IF ( first_time ) THEN
! compute the number of images for a given cut-off
ALLOCATE ( n_images ( nkind, nkind, 3 ), STAT = isos )
IF ( isos /= 0 ) CALL stop_memory ( 'overlap list', &
'n_images', 3 * nkind ** 2 )
CALL number_of_image_cells ( box, n_images, verlet_skin, rlist_cut )
! allocate the start of the lists
ALLOCATE (start_list_nl,STAT=isos)
IF ( isos /= 0 ) CALL stop_memory ( 'overlap_list', 'start_list', 0 )
NULLIFY (start_list_nl%next)
IF (MAXVAL(n_images) /= 0 ) THEN
npairs=natoms*(natoms+1)/2
ALLOCATE(start_list_im(npairs),STAT=isos)
IF ( isos /= 0 ) CALL stop_memory ( 'overlap_list', 'start_list_im', 0 )
DO ipair=1,npairs
ALLOCATE(start_list_im(ipair)%image, STAT=isos)
IF ( isos /= 0 ) CALL stop_memory ( 'overlap_list', &
'start_list_im%image', 0 )
NULLIFY (start_list_im(ipair)%image%next)
END DO
END IF
END IF
IF (first_time .AND. globenv % print_level > 4 &
.OR. globenv % print_level > 9 ) THEN
IF (globenv % ionode) WRITE (globenv % scr,'(A,T71,i10 )' ) &
' overlap_list| # of atoms : ', natoms
END IF
!
! let the current_neighbor point to the start of the list
!
current_neighbor => start_list_nl
!
! loop over all distinct i,j pairs of atoms
!
ipair=0
DO ikind = 1, nkind
nat = size (ki(ikind)%atom_list)
DO iatom = 1, nat
ii = ki(ikind)%atom_list(iatom)
i = part(ii)%node
pnode(i) %nsneighbor = 0
nexcl = pnode(i) %nexcl
IF (nexcl>size(elist)) STOP 'els'
ll_exclude => pnode(i) %ex
! list number of exclusion
DO iexclude = 1, nexcl
elist(iexclude) = ll_exclude%p%iatom
ll_exclude => ll_exclude%next
END DO
! sort elist
CALL sort(elist,nexcl,ind)
!
! pnode(i)%sl is the start of the OVERLAP neighbor link_list for
! atom i. Let it point to the current position in
! the link_list
!
pnode(i) %sl => current_neighbor
! first: loop over all particles in the same cell
DO jkind = 1, nkind
ll => part(ii) % cell % pTYPE ( jkind ) % ll
DO jpart = 1, part(ii) % cell % natoms(jkind )
j = ll%atom
! use black/white scheme to avoid double counting
IF (include_list(ii,j)) THEN
ipair=ipair+1
IF (MAXVAL(n_images(ikind,jkind,:)) /= 0 ) &
current_neighbor%image=>start_list_im(ipair)%image
! look for a match in the exclusion list (j in excl of pnode(i))
SELECT CASE (nexcl)
CASE DEFAULT
IF ( .NOT. get_match(j,elist,nexcl)) THEN
CALL update_verlet_list(n_images ( ikind, jkind, : ), &
part(j), j, &
pnode(i),hmat,h_inv,perd,rlist_cutsq(ikind,jkind ), &
current_neighbor,list_type)
END IF
CASE (0)
CALL update_verlet_list(n_images(ikind,jkind,:),part(j),j, &
pnode(i),hmat,h_inv,perd,rlist_cutsq(ikind,jkind ), &
current_neighbor,list_type)
CASE (1)
IF ( j/=elist(1)) THEN
CALL update_verlet_list(n_images(ikind,jkind,:), &
part(j),j, &
pnode(i),hmat,h_inv,perd,rlist_cutsq(ikind,jkind ), &
current_neighbor,list_type)
END IF
END SELECT
END IF
ll => ll%next
END DO
END DO
! second: loop over all neighbor cells
DO jkind = 1, nkind
IF (ALL(n_cell==1)) EXIT
ncells=SIZE (neighbor(ikind,jkind ) % index, 2)
DO k = 1, ncells
index(1) = part(ii) % cell_ol % neighbor(ikind,jkind ) %index(1,k)
index(2) = part(ii) % cell_ol % neighbor(ikind,jkind ) %index(2,k)
index(3) = part(ii) % cell_ol % neighbor(ikind,jkind ) %index(3,k)
ll => cell_ll(index(1),index(2),index(3)) % pTYPE ( jkind ) %ll
DO jpart = 1, cell_ll(index(1),index(2),index(3)) %natoms(jkind)
ipair=ipair+1
j = ll%atom
SELECT CASE (nexcl)
CASE DEFAULT
IF ( .NOT. get_match(j,elist,nexcl)) THEN
CALL update_verlet_list(n_images(ikind,jkind,:), &
part(j),j, &
pnode(i),hmat,h_inv,perd,rlist_cutsq(ikind,jkind ), &
current_neighbor,list_type)
END IF
CASE (0)
CALL update_verlet_list(n_images(ikind,jkind,:),part(j),j, &
pnode(i),hmat,h_inv,perd,rlist_cutsq(ikind,jkind ), &
current_neighbor,list_type)
CASE (1)
IF ( j/=elist(1)) THEN
CALL update_verlet_list(n_images(ikind,jkind,:), &
part(j),j, &
pnode(i),hmat,h_inv,perd,rlist_cutsq(ikind,jkind ), &
current_neighbor,list_type)
END IF
END SELECT
ll => ll%next
END DO
END DO
END DO
END DO
IF ( first_time .AND. globenv % print_level>4 &
.OR. globenv % print_level>9) THEN
IF (globenv % ionode) WRITE (globenv % scr,'(A,i10,T54,A,T71,I10 )' ) &
' overlap_list| Atom :', i, ' # of neighbors:', &
pnode(i) %nsneighbor
END IF
END DO
!
! write some info on the neighbor lists to output
!
nloc = sum(pnode ( : )%nsneighbor)
IF ( first_time .AND. globenv % print_level > 4 &
.OR. globenv % print_level > 9 ) THEN
#if defined(__parallel)
nsum = nloc
CALL mp_sum(nsum,globenv%group)
nmax = nloc
CALL mp_max(nmax,globenv%group)
nmin = nloc
CALL mp_min(nmin,globenv%group)
IF (globenv % ionode) WRITE (globenv % scr,'(A,T69,I12 )' ) &
' overlap_list| Total number of interactions ', nsum
IF (globenv % ionode) WRITE (globenv % scr,'(A,T69,I12 )' ) &
' overlap_list| Minimum number of interactions per node', nmin
IF (globenv % ionode) WRITE (globenv % scr,'(A,T69,I12 )' ) &
' overlap_list| Maximum number of interactions per node', nmax
IF (globenv % ionode) WRITE (globenv % scr,'(A,T69,F10.2,A )' ) &
' overlap_list| Load balancing error', 100.0_dbl &
* REAL ( nmax - nmin, dbl) &
/ REAL ( nsum, dbl ) * globenv % num_pe, ' %'
#else
IF (globenv % ionode) WRITE (globenv % scr,'(A,T69,I12 )' ) &
' overlap_list| Total number of interactions ', nloc
#endif
END IF
fis = nloc * 1.0E-3_dbl
CALL timestop(fis,handle)
DEALLOCATE (elist,STAT=isos)
IF ( isos /= 0 ) CALL stop_memory ( 'overlap_list', 'elist' )
DEALLOCATE (ind,STAT=isos)
IF ( isos /= 0 ) CALL stop_memory ( 'overlap_list', 'ind' )
END SUBROUTINE overlap_list
!******************************************************************************
END MODULE linklist_overlap_list

18
src/linklist_utilities.F
View file

@ -13,8 +13,8 @@
!! CJM
!!
!! MODIFICATION HISTORY
!! none
!!
!! gloria 8 may 2001 : linklist_verlet_list : added argument list_type,
!! adapt the routine to be used as well for the overlaps lists
!! SOURCE
!******************************************************************************
@ -63,7 +63,7 @@ MODULE linklist_utilities
END FUNCTION get_match
!------------------------------------------------------------------------------!
SUBROUTINE update_verlet_list(n_images,part,j,pnode,hmat,h_inv,perd, &
distsq,current_neighbor)
distsq,current_neighbor,list_type)
IMPLICIT NONE
TYPE (particle_type ), INTENT ( IN ), TARGET :: part
TYPE (particle_node_type ), INTENT ( INOUT ) :: pnode
@ -72,10 +72,13 @@ MODULE linklist_utilities
REAL ( dbl ), INTENT ( IN ) :: hmat(3,3), h_inv(3,3), perd(3)
REAL ( dbl ), INTENT ( IN ) :: distsq
TYPE (linklist_neighbor), POINTER :: current_neighbor
INTEGER, intent (in) :: list_type
! locals
REAL ( dbl ), DIMENSION (3) :: rij, s
REAL ( dbl ) :: rijsq
INTEGER :: isos
! initialize rijsq to large positive numer
rijsq = 1.e10_dbl
rij(1) = part % r(1) - pnode % p % r(1)
rij(2) = part % r(2) - pnode % p % r(2)
@ -105,7 +108,14 @@ MODULE linklist_utilities
! position in the link_list point to the jth atom
!
IF (rijsq<=distsq.OR.MAXVAL(n_images) /= 0 ) THEN
pnode%nneighbor = pnode%nneighbor + 1
SELECT CASE (list_type)
CASE (1)
pnode%nneighbor = pnode%nneighbor + 1
CASE (2)
pnode%nsneighbor = pnode%nsneighbor + 1
END SELECT
current_neighbor%p => part
current_neighbor%index = j
current_neighbor%nimages = 0

18
src/linklist_verlet_list.F
View file

@ -52,7 +52,7 @@ CONTAINS
SUBROUTINE verlet_list(natom_types,n_images,n_cell,pnode,part,box, &
verlet_skin,rlist_cut,rlist_cutsq,neighbor, &
cell_ll,start_list_nl,start_list_im, globenv)
cell_ll,start_list_nl,start_list_im, globenv,list_type)
IMPLICIT NONE
@ -67,9 +67,6 @@ SUBROUTINE verlet_list(natom_types,n_images,n_cell,pnode,part,box, &
TYPE ( cell_type ), INTENT ( IN ) :: box
TYPE ( linklist_neighbor), POINTER :: start_list_nl
TYPE ( startlist_images_type ), POINTER :: start_list_im ( : )
!*apsi* TYPE ( cell_neighborlist_type ), INTENT ( IN ) :: cell_ll(:,:,:)
!*apsi* TYPE ( neighbor_cells_type ), INTENT ( IN ) :: neighbor ( :, : )
!*apsi* PGI 3.2 did not like these as arrays
TYPE ( cell_neighborlist_type ), DIMENSION ( :, :, : ), POINTER :: cell_ll
TYPE ( neighbor_cells_type ), DIMENSION ( :, : ), POINTER :: neighbor
TYPE ( global_environment_type ), INTENT ( IN ) :: globenv
@ -78,6 +75,7 @@ SUBROUTINE verlet_list(natom_types,n_images,n_cell,pnode,part,box, &
INTEGER, DIMENSION ( :, :, : ), POINTER :: n_images
INTEGER, DIMENSION ( : ), INTENT ( IN ) :: n_cell
INTEGER, INTENT ( IN ) :: natom_types
INTEGER, INTENT (IN) :: list_type
! Locals
TYPE (cell_atoms_linklist_type ), POINTER :: ll
@ -213,18 +211,18 @@ SUBROUTINE verlet_list(natom_types,n_images,n_cell,pnode,part,box, &
CALL update_verlet_list(n_images ( itype, jtype, : ), &
part(j), j, &
pnode(i),hmat,h_inv,perd,rlist_cutsq(itype,jtype ), &
current_neighbor)
current_neighbor,list_type)
END IF
CASE (0)
CALL update_verlet_list(n_images(itype,jtype,:),part(j),j, &
pnode(i),hmat,h_inv,perd,rlist_cutsq(itype,jtype ), &
current_neighbor)
current_neighbor,list_type)
CASE (1)
IF ( j/=elist(1)) THEN
CALL update_verlet_list(n_images(itype,jtype,:), &
part(j),j, &
pnode(i),hmat,h_inv,perd,rlist_cutsq(itype,jtype ), &
current_neighbor)
current_neighbor,list_type)
END IF
END SELECT
END IF
@ -251,18 +249,18 @@ SUBROUTINE verlet_list(natom_types,n_images,n_cell,pnode,part,box, &
CALL update_verlet_list(n_images(itype,jtype,:), &
part(j),j, &
pnode(i),hmat,h_inv,perd,rlist_cutsq(itype,jtype ), &
current_neighbor)
current_neighbor,list_type)
END IF
CASE (0)
CALL update_verlet_list(n_images(itype,jtype,:),part(j),j, &
pnode(i),hmat,h_inv,perd,rlist_cutsq(itype,jtype ), &
current_neighbor)
current_neighbor,list_type)
CASE (1)
IF ( j/=elist(1)) THEN
CALL update_verlet_list(n_images(itype,jtype,:), &
part(j),j, &
pnode(i),hmat,h_inv,perd,rlist_cutsq(itype,jtype ), &
current_neighbor)
current_neighbor,list_type)
END IF
END SELECT
ll => ll%next

110
src/md.F
View file

@ -39,20 +39,21 @@ MODULE md
TYPE simulation_parameters_type
INTEGER :: nsteps
REAL ( dbl ) :: dt
REAL ( dbl ) :: temp_ext
REAL ( dbl ) :: temp_ext, pol_temp
REAL ( dbl ) :: p_ext
CHARACTER ( LEN = 20 ) :: ensemble
CHARACTER ( LEN = 20 ) :: read_type
CHARACTER ( LEN = 20 ) :: ensemble, pol_ensemble
CHARACTER ( LEN = 20 ) :: read_type, pol_read_type
CHARACTER ( LEN = 20 ) :: program
CHARACTER ( LEN = 20 ) :: pol_option
LOGICAL :: constraint
INTEGER :: nfree
INTEGER :: nc
INTEGER :: nyosh
INTEGER :: nhclen
INTEGER :: nfree, pol_nfree
INTEGER :: nc, pol_nc
INTEGER :: nyosh, pol_nyosh
INTEGER :: nhclen, pol_nhclen
INTEGER, DIMENSION ( 3 ) :: n_cell
REAL ( dbl ) :: tau_nhc
REAL ( dbl ) :: tau_nhc, pol_tau_nhc
REAL ( dbl ) :: tau_cell
REAL ( dbl ), POINTER, DIMENSION ( : ) :: dt_yosh
REAL ( dbl ), POINTER, DIMENSION ( : ) :: dt_yosh, pol_dt_yosh
REAL ( dbl ) :: shake_tol
REAL ( dbl ) :: verlet_skin
TYPE ( pimd_parameters_type ) :: pimd_params
@ -98,14 +99,22 @@ CONTAINS
!!I SECTION: &md ... &end I
!!I I
!!I ensemble: [nve,nvt] ensemble type I
!!I pol ensemble: [nve,nvt] ensemble type I
!!I steps: n total number of steps I
!!I timestep: dt time step [fs] I
!!I temperature: T temperature [K] I
!!I coef temperature: T temperature [K] I
!!I restart: [pos,posvel] start type I
!!I constraints: [on,off] constraints I
!!I shake: eps convergence for Shake algorithm I
!!I verlet_skin vs buffer size for Verlet lists I
!!I cells x y z Number of cells in each direction I
!!I Pol_Nose_parameter & parameter of Nose-Hoover-chains I
!!I option: pol_option & thermostat option I
!!I length: nhclen & length of thermostat I
!!I Yoshida: nyosh & order of Yoshida integrator I
!!I mts: nc & number of multiple time steps I
!!I timecon: tau_nhc NHC time constant I
!!I Nose_parameter & parameter of Nose-Hoover-chains I
!!I length: nhclen & length of thermostat I
!!I Yoshida: nyosh & order of Yoshida integrator I
@ -155,17 +164,25 @@ SUBROUTINE read_md_section ( simpar, mdpar, mdio )
!..defaults
simpar % ensemble = 'NVE'
simpar % pol_ensemble = 'NVE'
simpar % pol_option = 'GLOBAL'
simpar % nsteps = 0
simpar % dt = 0.0_dbl
simpar % temp_ext = 0.0_dbl
simpar % pol_temp = 0.0_dbl
simpar % p_ext = 0.0_dbl
simpar % read_type = 'POS'
simpar % pol_read_type = 'INIT'
simpar % constraint = .FALSE.
simpar % shake_tol = 1.0E-6_dbl
simpar % nhclen = 1
simpar % pol_nhclen = 1
simpar % pol_nc = 1
simpar % nc = 1
simpar % n_cell = 1
simpar % nyosh = 1
simpar % pol_nyosh = 1
simpar % pol_tau_nhc = 1000.0_dbl
simpar % tau_nhc = 1000.0_dbl
simpar % tau_cell = 1000.0_dbl
simpar % verlet_skin = 1.0_dbl
@ -214,6 +231,12 @@ SUBROUTINE read_md_section ( simpar, mdpar, mdio )
CALL uppercase ( string )
simpar % ensemble = string ( 1:ilen )
CASE ( 'POL_E' )
ilen = 0
CALL cfield ( string, ilen )
CALL uppercase ( string )
simpar % pol_ensemble = string ( 1:ilen )
CASE ( 'STEPS' )
simpar % nsteps = get_int()
@ -231,12 +254,21 @@ SUBROUTINE read_md_section ( simpar, mdpar, mdio )
CASE ( 'TEMPER' )
simpar % temp_ext = get_real()
CASE ( 'POL_T' )
simpar % pol_temp = get_real()
CASE ( 'RESTAR' )
ilen = 0
CALL cfield ( string, ilen )
CALL uppercase ( string )
simpar % read_type = string ( 1:ilen )
CASE ( 'POLRES' )
ilen = 0
CALL cfield ( string, ilen )
CALL uppercase ( string )
simpar % pol_read_type = string ( 1:ilen )
CASE ( 'CONSTR' )
ilen = 0
CALL cfield ( string, ilen )
@ -247,6 +279,34 @@ SUBROUTINE read_md_section ( simpar, mdpar, mdio )
CASE ( 'SHAKE' )
simpar % shake_tol = get_real()
CASE ( 'POL_N' )
DO WHILE ( test_next() == 'C' )
ilen = 6
CALL cfield ( string, ilen )
CALL uppercase ( string )
SELECT CASE ( string )
CASE DEFAULT
CALL p_error()
CALL stop_parser ( 'read_md_section%nose', &
'unknown suboption' )
CASE ( 'OPTION' )
ilen = 0
CALL cfield ( string, ilen )
CALL uppercase ( string )
simpar % pol_option = string ( 1:ilen )
CASE ( 'LENGTH' )
simpar % pol_nhclen = get_int()
CASE ( 'YOSHID' )
simpar % pol_nyosh = get_int()
CASE ( 'TIMECO' )
simpar % pol_tau_nhc = get_real()
CASE ( 'MTS' )
simpar % pol_nc = get_int()
END SELECT
END DO
CASE ( 'NOSE_P' )
DO WHILE ( test_next() == 'C' )
@ -414,6 +474,21 @@ SUBROUTINE read_md_section ( simpar, mdpar, mdio )
WRITE ( iw, '( A,T71,G10.4 )' ) ' MD| Tolerance for shake ', &
simpar % shake_tol
END IF
IF ( simpar % pol_ensemble /= 'NVE' ) THEN
WRITE ( iw, '( A,T61,A )' ) ' MD| Pol Ensemble Type ', &
ADJUSTR ( simpar % pol_ensemble )
WRITE ( iw, '( A,T71,F10.2 )' ) ' MD| Pol Temperature [K] ', &
simpar % pol_temp
WRITE ( iw, '( A )' ) ' MD| Pol Nose-Hoover-Chain parameters '
WRITE ( iw, '( A,T77,I4 )' ) ' MD| Pol Nose-Hoover-Chain length ', &
simpar % pol_nhclen
WRITE ( iw, '( A,T71,F10.4 )' ) &
' MD| Pol Nose-Hoover-Chain time constant [fs] ', simpar % pol_tau_nhc
WRITE ( iw, '( A,T77,I4 )' ) ' MD| Pol Order of Yoshida integrator ', &
simpar % pol_nyosh
WRITE ( iw, '( A,T77,I4 )' ) ' MD| Pol Number of multiple time steps', &
simpar % pol_nc
ENDIF
WRITE ( iw, '( A,T63,i10,A )' ) ' MD| Print MD information every', &
mdio % iscreen, ' step(s)'
WRITE ( iw, '( A,T20,A,T71,A )' ) ' MD| File type', &
@ -488,16 +563,19 @@ SUBROUTINE initialize_velocities ( simpar, part, mdpar )
! initializing velocities
IF (mdpar % ionode) THEN
DO i = 1, natoms
part(i) % v(1) = gasdev(mdpar % idum)
part(i) % v(2) = gasdev(mdpar % idum)
part(i) % v(3) = gasdev(mdpar % idum)
part ( i ) % v ( 1 ) = gasdev ( mdpar % idum )
part ( i ) % v ( 2 ) = gasdev ( mdpar % idum )
part ( i ) % v ( 3 ) = gasdev ( mdpar % idum )
END DO
akin = 0.0_dbl
DO i = 1, natoms
mass = part(i) % prop % mass
akin = akin + 0.5_dbl*(mass*part(i) % v(1)*part(i) % v(1))
akin = akin + 0.5_dbl*(mass*part(i) % v(2)*part(i) % v(2))
akin = akin + 0.5_dbl*(mass*part(i) % v(3)*part(i) % v(3))
mass = part ( i ) % prop % mass
akin = akin + 0.5_dbl * &
( mass * part ( i ) % v ( 1 ) * part ( i ) % v ( 1 ) )
akin = akin + 0.5_dbl * &
( mass * part ( i ) % v ( 2 ) * part ( i ) % v ( 2 ) )
akin = akin + 0.5_dbl * &
( mass * part ( i ) % v ( 3 ) * part ( i ) % v ( 3 ) )
END DO
temp = 2.0_dbl * akin / REAL ( simpar % nfree, dbl )

8
src/method_specifications.F
View file

@ -38,7 +38,7 @@ MODULE method_specifications
! *** Public subroutines ***
PUBLIC :: read_method_specifications,&
write_method_specifications
write_method_specifications, get_maxder
! *****************************************************************************
@ -179,6 +179,12 @@ CONTAINS
END SUBROUTINE write_method_specifications
! *****************************************************************************
SUBROUTINE get_maxder(maxder_index)
IMPLICIT NONE
INTEGER, INTENT(in), OPTIONAL :: maxder_index ! maxder
IF (PRESENT(maxder_index)) maxder = maxder_index
END SUBROUTINE get_maxder
! *****************************************************************************
END MODULE method_specifications

74
src/molecule_input.F
View file

@ -235,7 +235,7 @@ SUBROUTINE read_molecule_section ( mol_setup, set_fn, globenv )
! Locals
INTEGER :: ierror, ilen, iw, source, group, icount, nmol_type, ios
INTEGER :: i, j, imol, na, nat, j1, j2, nf, nt, ie, ia, i1, i2, i3, i4
INTEGER :: isos, nbasis, ndrho, allgrp
INTEGER :: isos, nbasis, ndrho, allgrp, npol, nrho0
CHARACTER ( LEN = 20 ) :: libstring, string, string2, at1, at2, at3, at4
CHARACTER ( LEN = 6 ) :: label
CHARACTER ( LEN = 20 ), DIMENSION ( : ), ALLOCATABLE :: molname
@ -277,10 +277,12 @@ SUBROUTINE read_molecule_section ( mol_setup, set_fn, globenv )
mol_setup(imol) %molpar%ntorsions = 0
mol_setup(imol) %molpar%ndcon = 0
mol_setup(imol) %molpar%n3x3con = 0
mol_setup(imol) %molpar%nbasis_pol = 0
mol_setup(imol) %molpar%nbasis_rho0 = 0
DO
CALL read_line
ilen = 6
ilen = 8
CALL cfield(string2,ilen)
CALL uppercase ( string2 )
@ -387,6 +389,39 @@ SUBROUTINE read_molecule_section ( mol_setup, set_fn, globenv )
mol_setup(imol) %molpar%n3x3c(3,i) = get_int()
END DO
END IF
CASE ('RHO0_SET')
na = mol_setup(imol) %molpar%natom
IF (.NOT. ASSOCIATED (mol_setup(imol)%molpar%rho0_basis)) THEN
ALLOCATE (mol_setup(imol)%molpar%rho0_basis(na),STAT=ios)
IF ( ios /= 0 ) CALL stop_memory &
( 'read_molecule_section', 'rho0_basis', na )
END IF
nrho0 = mol_setup(imol)%molpar%nbasis_rho0 + 1
mol_setup(imol)%molpar%nbasis_rho0 = nrho0
ilen=0
CALL cfield(mol_setup(imol)%molpar%rho0_basis(nrho0)%aname,ilen)
ilen=0
CALL cfield(mol_setup(imol)%molpar%rho0_basis(nrho0)%libname,ilen)
mol_setup(imol)%molpar%rho0_basis(nrho0)%nmax=get_int()
mol_setup(imol)%molpar%rho0_basis(nrho0)%eps=get_real()
CASE ('POL_SET')
na = mol_setup(imol) %molpar%natom
IF (.NOT. ASSOCIATED (mol_setup(imol)%molpar%pol_basis)) THEN
ALLOCATE (mol_setup(imol)%molpar%pol_basis(na),STAT=ios)
IF ( ios /= 0 ) CALL stop_memory &
( 'read_molecule_section', 'pol_basis', na )
END IF
npol = mol_setup(imol)%molpar%nbasis_pol + 1
mol_setup(imol)%molpar%nbasis_pol = npol
ilen=0
CALL cfield(mol_setup(imol)%molpar%pol_basis(npol)%aname,ilen)
ilen=0
CALL cfield(mol_setup(imol)%molpar%pol_basis(npol)%libname,ilen)
mol_setup(imol)%molpar%pol_basis(npol)%nmax=get_int()
mol_setup(imol)%molpar%pol_basis(npol)%eps=get_real()
CASE ( 'END')
ilen = 0
@ -554,6 +589,41 @@ SUBROUTINE read_molecule_section ( mol_setup, set_fn, globenv )
IF ( isos /= 0 ) CALL stop_memory ( 'read_molecule_section', &
'tuple' )
END IF
na = mol_setup(imol) %molpar%nbasis_rho0
IF (na>0) THEN
WRITE ( iw, '( A,T71,I10 )' ) &
' MOLECULE| Number of effective densities', &
mol_setup(imol) %molpar%nbasis_rho0
DO i = 1, mol_setup(imol)%molpar%nbasis_rho0
WRITE (iw,'(A,T61,A)') ' MOLECULE| name of rho0 basis set ', &
mol_setup(imol)%molpar%rho0_basis(i)%aname
WRITE (iw,'(A,T61,A)') ' MOLECULE| library for basis set ', &
mol_setup(imol)%molpar%rho0_basis(i)%libname
WRITE (iw,'(A,T71,I10)') ' MOLECULE| max grid points ', &
mol_setup(imol)%molpar%rho0_basis(i)%nmax
WRITE (iw,'(A,T68,F12.6)') ' MOLECULE| epsilon ', &
mol_setup(imol)%molpar%rho0_basis(i)%eps
END DO
END IF
na = mol_setup(imol) %molpar%nbasis_pol
IF (na>0) THEN
WRITE ( iw, '( A,T71,I10 )' ) &
' MOLECULE| Number of polarization densities', &
mol_setup(imol) %molpar%nbasis_pol
DO i = 1, mol_setup(imol)%molpar%nbasis_pol
WRITE (iw,'(A,T61,A)') ' MOLECULE| name of polarizable basis set ', &
mol_setup(imol)%molpar%pol_basis(i)%aname
WRITE (iw,'(A,T61,A)') ' MOLECULE| library for basis set ', &
mol_setup(imol)%molpar%pol_basis(i)%libname
WRITE (iw,'(A,T71,I10)') ' MOLECULE| max grid points ', &
mol_setup(imol)%molpar%pol_basis(i)%nmax
WRITE (iw,'(A,T68,F12.6)') ' MOLECULE| epsilon ', &
mol_setup(imol)%molpar%pol_basis(i)%eps
END DO
END IF
WRITE ( iw,'( )' )
END DO
END IF

13
src/molecule_types.F
View file

@ -32,7 +32,7 @@ MODULE molecule_types
linklist_g3x3_constraints, linklist_neighbor, linklist_exclusion, &
particle_node_type, localization_list, linklist_images
PUBLIC :: particle_type, particle_prop_type
PUBLIC :: molecule_type, molecule_parameter_type
PUBLIC :: molecule_type, molecule_parameter_type, read_basis_type
PUBLIC :: intra_parameters_type
TYPE linklist_images
@ -137,6 +137,13 @@ MODULE molecule_types
REAL ( dbl ), POINTER :: constraint_distance ( :, : )
END TYPE intra_parameters_type
TYPE read_basis_type
CHARACTER (len=20) :: aname
CHARACTER (len=20) :: libname
REAL(dbl) :: eps
INTEGER :: nmax
END TYPE read_basis_type
TYPE molecule_parameter_type
CHARACTER ( LEN = 20 ) :: name
INTEGER :: natom
@ -155,6 +162,10 @@ MODULE molecule_types
INTEGER, POINTER :: ndc ( :, : )
INTEGER :: n3x3con
INTEGER, POINTER :: n3x3c ( :, : )
INTEGER :: nbasis_pol
TYPE(read_basis_type), POINTER :: pol_basis(:)
INTEGER :: nbasis_rho0
TYPE(read_basis_type), POINTER :: rho0_basis(:)
END TYPE molecule_parameter_type
TYPE molecule_type

215
src/nhc.F
View file

@ -15,12 +15,14 @@
!! MODIFICATION HISTORY
!! CJM, 20-Feb-01
!! JGH (10-Mar-2001)
!! CJM (10-Apr-2001)
!!
!! SOURCE
!******************************************************************************
MODULE nhc
USE coefficient_types, ONLY : coeff_type
USE global_types, ONLY : global_environment_type
USE kinds, ONLY : dbl
USE md, ONLY : simulation_parameters_type
@ -33,16 +35,123 @@ MODULE nhc
IMPLICIT NONE
PRIVATE
PUBLIC :: nhc_control_particles, nhc_control_barostat
REAL ( dbl ), TARGET :: scale
PUBLIC :: nhc_particle_setup, nhc_barostat_setup, nhc_pol_setup
!!*****
!******************************************************************************
CONTAINS
!******************************************************************************
!!****** nhc/nhc_control_barostat [1.0] *
!!****** nhc/nhc_pol_setup [1.0] *
!!
!! NAME
!! nhc_pol_setup
!!
!! FUNCTION
!! Creates the thermostatting for polarization degrees of freedom
!!
!! AUTHOR
!! CJM 10-Apr-01
!!
!! MODIFICATION HISTORY
!! NONE
!!
!! SOURCE
!******************************************************************************
SUBROUTINE nhc_pol_setup ( simpar, coeff, nhc, globenv )
! creates the thermostatting maps use global_mod
IMPLICIT NONE
! begin global declaratation
TYPE ( simulation_parameters_type ), INTENT ( INOUT ) :: simpar
TYPE ( coeff_type ), INTENT ( IN ) :: coeff
TYPE ( lnhc_parameters_type ), INTENT ( OUT ) :: nhc
TYPE ( global_environment_type ), INTENT ( IN ) :: globenv
! begin local definitions
INTEGER :: i, iw
INTEGER :: isos, number
!
iw = globenv % scr
!
SELECT CASE ( simpar % pol_ensemble ( 1 : 3 ) )
CASE DEFAULT
ALLOCATE ( nhc % nvt ( 0, 0 ), STAT = isos )
IF (isos/=0) CALL stop_program('nhc_ao_setup', &
'failed to allocate nhc%nvt' )
nhc % nyosh = 0
nhc % nc = 0
nhc % nhc_len = 0
nhc % num_nhc = 0
CASE ( 'NVT', 'NPT' )
nhc % nyosh = simpar % pol_nyosh
nhc % nc = simpar % pol_nc
nhc % nhc_len = simpar % pol_nhclen
IF ( simpar % pol_option == 'GLOBAL' ) THEN
number = 1
nhc % num_nhc = 1
ELSEIF ( simpar % pol_option == 'MASSIVE' ) THEN
number = SIZE ( coeff % ao % cr )
nhc % num_nhc = number
ENDIF
! Now we can allocate the target array s_kin and p_kin
ALLOCATE ( nhc % s_kin ( nhc % num_nhc ), STAT = isos)
IF ( isos /= 0 ) CALL stop_program ( 'nhc_ao_setup', &
'failed to allocate nhc%s_kin' )
ALLOCATE ( nhc % v_scale ( nhc % num_nhc ), STAT = isos )
IF (isos/=0) CALL stop_program('nhc_ao_setup', &
'failed to allocate nhc%v_scale' )
ALLOCATE ( nhc % p_kin ( 1, number ), STAT = isos )
IF (isos/=0) CALL stop_program('nhc_ao_setup', &
'failed to allocate nhc%p_kin' )
ALLOCATE ( nhc % p_scale ( 1, number ), STAT = isos )
IF (isos/=0) CALL stop_program('nhc_ao_setup', &
'failed to allocate nhc%p_scale')
! Begin the mapping loop
IF ( simpar % pol_option == 'GLOBAL' ) THEN
DO i = 1, number
nhc % p_kin ( 1, i ) % point => nhc % s_kin ( 1 )
nhc % p_scale ( 1, i ) % point => nhc % v_scale ( 1 )
END DO
ELSEIF ( simpar % pol_option == 'MASSIVE' ) THEN
DO i = 1, number
nhc % p_kin ( 1, i ) % point => nhc % s_kin ( i )
nhc % p_scale ( 1, i ) % point => nhc % v_scale ( i )
END DO
ENDIF
ALLOCATE ( nhc % nvt ( nhc % nhc_len, nhc % num_nhc ), STAT = isos )
IF (isos/=0) CALL stop_program('nhc_ao_setup', &
'failed to allocate nhc%nvt' )
!
! Now that we know how many there are stick this into nhc % static % nkt
! (number of degrees of freedom times k_B T )
DO i = 1, nhc % num_nhc
nhc % nvt ( 1, i ) % nkt = simpar % pol_temp * number
END DO
! getting the number of degrees of freedom times k_B T for the rest of the chain
DO i = 2, nhc % nhc_len
nhc % nvt ( i, : ) % nkt = simpar % pol_temp
END DO
! Set up the Yoshida weights
ALLOCATE (nhc % dt_yosh (1 : nhc % nyosh ), STAT = isos )
IF (isos/=0) CALL stop_program('nhc_pol_setup', &
'failed to allocate nhc%dt_yosh')
CALL yoshida_coef( nhc, simpar % dt )
END SELECT
END SUBROUTINE nhc_pol_setup
!******************************************************************************
!!****** nhc/nhc_barostat_setup [1.0] *
!!
!! NAME
!! nhc_barostat
@ -60,7 +169,7 @@ CONTAINS
!! SOURCE
!******************************************************************************
SUBROUTINE nhc_control_barostat ( simpar, nhc, globenv )
SUBROUTINE nhc_barostat_setup ( simpar, nhc, globenv )
! creates the thermostatting maps use global_mod
IMPLICIT NONE
! begin global declaratation
@ -69,14 +178,14 @@ SUBROUTINE nhc_control_barostat ( simpar, nhc, globenv )
TYPE ( global_environment_type ), INTENT ( IN ) :: globenv
! begin local definitions
INTEGER :: i, iw
INTEGER :: isos, ncoef
INTEGER :: isos, number
!
iw = globenv % scr
!
SELECT CASE ( simpar % ensemble ( 1 : 3 ) )
CASE DEFAULT
ALLOCATE ( nhc % nvt ( 0, 0 ), STAT = isos )
IF (isos/=0) CALL stop_program('nhc_control', &
IF (isos/=0) CALL stop_program('nhc_barostat_setup', &
'failed to allocate nhc%nvt' )
nhc % nyosh = 0
nhc % nc = 0
@ -93,44 +202,42 @@ SUBROUTINE nhc_control_barostat ( simpar, nhc, globenv )
nhc % num_nhc = 1
IF ( simpar % ensemble == 'NPT_F' ) THEN
ncoef = 9
number = 9
ELSE
ncoef = 1
number = 1
ENDIF
! Now we can allocate the target array s_kin and p_kin
ALLOCATE ( nhc % s_kin ( nhc % num_nhc ), STAT = isos)
IF ( isos /= 0 ) CALL stop_program ( 'nhc_control', &
IF ( isos /= 0 ) CALL stop_program ( 'nhc_barostat_setup', &
'failed to allocate nhc%s_kin' )
ALLOCATE ( nhc % v_scale ( nhc % num_nhc ), STAT = isos )
IF (isos/=0) CALL stop_program('nhc_control', &
IF (isos/=0) CALL stop_program('nhc_barostat_setup', &
'failed to allocate nhc%v_scale' )
ALLOCATE ( nhc % p_kin ( 1, ncoef ), STAT = isos )
IF (isos/=0) CALL stop_program('nhc_control', &
ALLOCATE ( nhc % p_kin ( 1, number ), STAT = isos )
IF (isos/=0) CALL stop_program('nhc_barostat_setup', &
'failed to allocate nhc%p_kin' )
ALLOCATE ( nhc % p_scale ( 1, ncoef ), STAT = isos )
IF (isos/=0) CALL stop_program('nhc_control', &
ALLOCATE ( nhc % p_scale ( 1, number ), STAT = isos )
IF (isos/=0) CALL stop_program('nhc_barostat_setup', &
'failed to allocate nhc%p_scale')
! initialize variable 'scale' for the 'none' option
scale = 1._dbl
! Begin the mapping loop
DO i = 1, ncoef
DO i = 1, number
nhc % p_kin ( 1, i ) % point => nhc % s_kin ( 1 )
nhc % p_scale ( 1, i ) % point => nhc % v_scale ( 1 )
END DO
ALLOCATE ( nhc % nvt ( nhc % nhc_len, nhc % num_nhc ), STAT = isos )
IF (isos/=0) CALL stop_program('nhc_control', &
IF (isos/=0) CALL stop_program('nhc_barostat_setup', &
'failed to allocate nhc%nvt' )
!
! Now that we know how many there are stick this into nhc % static % nkt
! (number of degrees of freedom times k_B T )
DO i = 1, nhc % num_nhc
nhc % nvt ( 1, i ) % nkt = simpar % temp_ext * ncoef
nhc % nvt ( 1, i ) % nkt = simpar % temp_ext * number
END DO
! getting the number of degrees of freedom times k_B T for the rest of the chain
@ -142,18 +249,18 @@ SUBROUTINE nhc_control_barostat ( simpar, nhc, globenv )
! Set up the Yoshida weights
ALLOCATE (nhc % dt_yosh (1 : nhc % nyosh ), STAT = isos )
IF (isos/=0) CALL stop_program('nhc_control', &
IF (isos/=0) CALL stop_program('nhc_barostat_setup', &
'failed to allocate nhc%dt_yosh')
CALL yoshida_coef( nhc, simpar % dt )
END SELECT
END SUBROUTINE nhc_control_barostat
END SUBROUTINE nhc_barostat_setup
!******************************************************************************
!!****** nhc/nhc_control_particles [1.0] *
!!****** nhc/nhc_particle_setup [1.0] *
!!
!! NAME
!! nhc_control_particles
!! nhc_particle_setup
!!
!! FUNCTION
!! Creates the thermostatting maps
@ -169,7 +276,7 @@ END SUBROUTINE nhc_control_barostat
!! SOURCE
!******************************************************************************
SUBROUTINE nhc_control_particles ( simpar, molecule, mol_setup, nhc, globenv )
SUBROUTINE nhc_particle_setup ( simpar, molecule, mol_setup, nhc, globenv )
IMPLICIT NONE
@ -191,6 +298,7 @@ SUBROUTINE nhc_control_particles ( simpar, molecule, mol_setup, nhc, globenv )
INTEGER, DIMENSION ( : ), ALLOCATABLE :: const_mol, tot_const
INTEGER, DIMENSION ( : ), ALLOCATABLE :: natoms_mol, mol_off
INTEGER, DIMENSION ( :, : ), ALLOCATABLE :: point
REAL ( dbl ), SAVE, TARGET :: scale
LOGICAL :: match, noneflag, nointer
!------------------------------------------------------------------------------
@ -200,7 +308,7 @@ SUBROUTINE nhc_control_particles ( simpar, molecule, mol_setup, nhc, globenv )
SELECT CASE ( simpar % ensemble ( 1 : 3 ) )
CASE ( 'NVE' )
ALLOCATE ( nhc % nvt ( 0, 0 ), STAT = isos )
IF ( isos /= 0 ) CALL stop_memory ( 'nhc_control', 'nhc%nvt', 0 )
IF ( isos /= 0 ) CALL stop_memory ( 'nhc_particle_setup', 'nhc%nvt', 0 )
nhc % nyosh = 0
nhc % nc = 0
nhc % nhc_len = 0
@ -210,23 +318,23 @@ SUBROUTINE nhc_control_particles ( simpar, molecule, mol_setup, nhc, globenv )
nhc%nc = simpar%nc
nhc%nhc_len = simpar%nhclen
! allocating the temporary arrays mnum and tnum, and nfree
! allocating the temporary arrays mnum and tnum
nmol_type = size(mol_setup)
nmol = size(molecule)
ALLOCATE (mnum(nmol_type ),STAT=isos)
IF ( isos /= 0 ) CALL stop_memory ( 'nhc_control', 'mnum', nmol_type )
IF ( isos /= 0 ) CALL stop_memory ( 'nhc_particle_setup', 'mnum', nmol_type )
ALLOCATE (tnum(nmol_type ),STAT=isos)
IF ( isos /= 0 ) CALL stop_memory ( 'nhc_control', 'tnum', nmol_type )
IF ( isos /= 0 ) CALL stop_memory ( 'nhc_particle_setup', 'tnum', nmol_type )
ALLOCATE (const_mol(nmol),STAT=isos)
IF ( isos /= 0 ) CALL stop_memory ( 'nhc_control', 'const_mol', nmol )
IF ( isos /= 0 ) CALL stop_memory ( 'nhc_particle_setup', 'const_mol', nmol )
ALLOCATE (tot_const(nmol),STAT=isos)
IF ( isos /= 0 ) CALL stop_memory ( 'nhc_control', 'tot_const', nmol )
IF ( isos /= 0 ) CALL stop_memory ( 'nhc_particle_setup', 'tot_const', nmol )
ALLOCATE (natoms_mol(nmol),STAT=isos)
IF ( isos /= 0 ) CALL stop_memory ( 'nhc_control', 'natoms_mol', nmol )
IF ( isos /= 0 ) CALL stop_memory ( 'nhc_particle_setup', 'natoms_mol', nmol )
ALLOCATE (point(2,nmol),STAT=isos)
IF ( isos /= 0 ) CALL stop_memory ( 'nhc_control','point', 2 * nmol )
IF ( isos /= 0 ) CALL stop_memory ( 'nhc_particle_setup','point', 2 * nmol )
ALLOCATE (mol_off(nmol_type+1),STAT=isos)
IF ( isos /= 0 ) CALL stop_memory ( 'nhc_control','mol_off', nmol_type + 1 )
IF ( isos /= 0 ) CALL stop_memory ( 'nhc_particle_setup','mol_off', nmol_type + 1 )
! Counting the total number of thermostats
number = 0
@ -279,7 +387,7 @@ SUBROUTINE nhc_control_particles ( simpar, molecule, mol_setup, nhc, globenv )
END IF
END DO
IF ( number == 0 ) THEN
CALL stop_program ( 'nhc_control', "Number of NHC = 0 -> use NVE" )
CALL stop_program ( 'nhc_paticle_setup', "Number of NHC = 0 -> use NVE" )
END IF
!
@ -304,7 +412,7 @@ SUBROUTINE nhc_control_particles ( simpar, molecule, mol_setup, nhc, globenv )
ELSE IF (mol_setup(j)%nhcopt=="MASSIVE") THEN
number = number + 3*mol_setup(j) %molpar%natom
ELSE
CALL stop_program ( 'nhc_control', '' )
CALL stop_program ( 'nhc_particle_setup', '' )
END IF
END DO
ELSE
@ -365,22 +473,22 @@ SUBROUTINE nhc_control_particles ( simpar, molecule, mol_setup, nhc, globenv )
SELECT CASE (simpar%ensemble(1:3))
CASE ( 'NPT' )
ALLOCATE (deg_of_freedom(number+1),STAT=isos)
IF ( isos /= 0 ) CALL stop_memory ( 'nhc_control', &
IF ( isos /= 0 ) CALL stop_memory ( 'nhc_particle_setup', &
'deg_of_freedom', number + 1 )
CASE DEFAULT
ALLOCATE (deg_of_freedom(number),STAT=isos)
IF ( isos /= 0 ) CALL stop_memory ( 'nhc_control', &
IF ( isos /= 0 ) CALL stop_memory ( 'nhc_particle_setup', &
'deg_of_freedom', number )
END SELECT
ALLOCATE (nhc%s_kin(number),STAT=isos)
IF ( isos /= 0 ) CALL stop_memory ( 'nhc_control', 'nhc%s_kin', number )
IF ( isos /= 0 ) CALL stop_memory ( 'nhc_particle_setup', 'nhc%s_kin', number )
ALLOCATE (nhc%v_scale(number),STAT=isos)
IF ( isos /= 0 ) CALL stop_memory ( 'nhc_control', 'nhc%v_scale', number )
IF ( isos /= 0 ) CALL stop_memory ( 'nhc_particle_setup', 'nhc%v_scale', number )
ALLOCATE (nhc%p_kin(3,natoms),STAT=isos)
IF ( isos /= 0 ) CALL stop_memory ( 'nhc_control', 'nhc%p_kin', 3 * natoms )
IF ( isos /= 0 ) CALL stop_memory ( 'nhc_particle_setup', 'nhc%p_kin', 3 * natoms )
ALLOCATE (nhc%p_scale(3,natoms),STAT=isos)
IF ( isos /= 0 ) CALL stop_memory ( 'nhc_control', &
IF ( isos /= 0 ) CALL stop_memory ( 'nhc_particle_setup', &
'nhc%p_scale', 3 * natoms )
! Check to see if there is at least one global thermostatting option defined
@ -518,7 +626,7 @@ SUBROUTINE nhc_control_particles ( simpar, molecule, mol_setup, nhc, globenv )
! doing massive thermostat option
ELSE IF ( mol_setup ( j ) % nhcopt == "MASSIVE" ) THEN
IF ( tot_const ( i ) > 0 ) THEN
CALL stop_program ( 'nhc_control', &
CALL stop_program ( 'nhc_particle_setup', &
"Massive thermostats with constraints not possible" )
END IF
DO ii = point(1,i), point(2,i)
@ -666,7 +774,7 @@ SUBROUTINE nhc_control_particles ( simpar, molecule, mol_setup, nhc, globenv )
! doing massive thermostat option
ELSE IF ( mol_setup ( i ) % nhcopt == "MASSIVE" ) THEN
IF ( tot_const ( i ) > 0 ) THEN
CALL stop_program ( 'nhc_control', &
CALL stop_program ( 'nhc_particle_setup', &
"Massive thermostats with constraints not possible" )
END IF
DO ii = point(1,i), point(2,i)
@ -694,20 +802,20 @@ SUBROUTINE nhc_control_particles ( simpar, molecule, mol_setup, nhc, globenv )
! deallocating the temporary arrays mnum and tnum
DEALLOCATE (mnum,STAT=isos)
IF ( isos /= 0 ) CALL stop_memory ( 'nhc_control','mnum' )
IF ( isos /= 0 ) CALL stop_memory ( 'nhc_particle_setup','mnum' )
DEALLOCATE (tnum,STAT=isos)
IF ( isos /= 0 ) CALL stop_memory ( 'nhc_control','tnum' )
IF ( isos /= 0 ) CALL stop_memory ( 'nhc_particle_setup','tnum' )
DEALLOCATE (const_mol,STAT=isos)
IF ( isos /= 0 ) CALL stop_memory ( 'nhc_control', 'const_mol' )
IF ( isos /= 0 ) CALL stop_memory ( 'nhc_particle_setup', 'const_mol' )
DEALLOCATE (tot_const,STAT=isos)
IF ( isos /= 0 ) CALL stop_memory ( 'nhc_control', 'tot_const' )
IF ( isos /= 0 ) CALL stop_memory ( 'nhc_particle_setup', 'tot_const' )
DEALLOCATE (natoms_mol,STAT=isos)
IF ( isos /= 0 ) CALL stop_memory ( 'nhc_control', 'natoms_mol' )
IF ( isos /= 0 ) CALL stop_memory ( 'nhc_particle_setup', 'natoms_mol' )
DEALLOCATE (mol_off,STAT=isos)
IF ( isos /= 0 ) CALL stop_memory ( 'nhc_control', 'mol_off' )
IF ( isos /= 0 ) CALL stop_memory ( 'nhc_particle_setup', 'mol_off' )
ALLOCATE (nhc%nvt(nhc%nhc_len,number),STAT=isos)
IF ( isos /= 0 ) CALL stop_memory ( 'nhc_control','nhc%nvt', 0 )
IF ( isos /= 0 ) CALL stop_memory ( 'nhc_particle_setup','nhc%nvt', 0 )
! Sum up the number of degrees of freedom on each thermostat.
! first: initialize the target
@ -744,6 +852,7 @@ SUBROUTINE nhc_control_particles ( simpar, molecule, mol_setup, nhc, globenv )
! re-initializing simpar%nfree to zero because of multiple thermostats
simpar % nfree = 0._dbl
DO i = 1, number
nhc % nvt ( 1, i ) % nkt = simpar % temp_ext * ( nhc % s_kin ( i ) - &
deg_of_freedom ( i ) )
@ -759,14 +868,14 @@ SUBROUTINE nhc_control_particles ( simpar, molecule, mol_setup, nhc, globenv )
! Set up the Yoshida weights
ALLOCATE ( nhc % dt_yosh ( 1 : nhc%nyosh ), STAT = isos )
IF ( isos /= 0 ) CALL stop_memory ( 'nhc_control', &
IF ( isos /= 0 ) CALL stop_memory ( 'nhc_particle_setup', &
'nhc%dt_yosh', nhc % nyosh )
CALL yoshida_coef ( nhc, simpar % dt )
END SELECT
END SUBROUTINE nhc_control_particles
END SUBROUTINE nhc_particle_setup
!!*****
!******************************************************************************

151
src/nose.F
View file

@ -13,14 +13,15 @@
!! CJM
!!
!! MODIFICATION HISTORY
!! CJM 20-Feb-2001: Now npt_ifo is allocated to zero when
!! not used
!! CJM 20-Feb-2001: Now npt_ifo is allocated to zero when not used
!! CJM 11-apr-2001: adding routines to thermostat ao_type
!!
!! SOURCE
!******************************************************************************
MODULE nose
USE ao_types, ONLY : ao_type
USE global_types, ONLY : global_environment_type
USE kinds, ONLY : dbl
USE molecule_types, ONLY : particle_node_type
@ -32,7 +33,7 @@ MODULE nose
IMPLICIT NONE
PRIVATE
PUBLIC :: lnhc, initial_npt, initial_nhc
PUBLIC :: lnhc, initialize_npt, initialize_nhc
PUBLIC :: lnhc_parameters_type, nhc_info_type, npt_info_type
PUBLIC :: yoshida_coef
PUBLIC :: particle_kinetic_type, particle_scale_type
@ -72,13 +73,137 @@ MODULE nose
END TYPE lnhc_parameters_type
INTERFACE lnhc
MODULE PROCEDURE lnhc_particles, lnhc_barostat
MODULE PROCEDURE lnhc_particles, lnhc_barostat, lnhc_ao
END INTERFACE
!!*****
!******************************************************************************
CONTAINS
!******************************************************************************
!!****** nose/lnhc_ao [1.0] *
!!
!! NAME
!! lnhc_ao
!!
!! FUNCTION
!!
!! AUTHOR
!! CJM
!!
!! CREATION DATE
!! 11-Apr-2001
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
SUBROUTINE lnhc_ao ( nhc, ao, group )
IMPLICIT NONE
! Arguments
TYPE ( lnhc_parameters_type ), INTENT ( INOUT ) :: nhc
TYPE ( ao_type ), INTENT ( INOUT ) :: ao
INTEGER, INTENT ( IN ) :: group
! Locals
INTEGER :: i, j, nnodes, inc, iyosh, n, inhc, handle, nx1, nx2, ncoef
REAL ( dbl ) :: scale, flops
CALL timeset('LNHC','I','Mflops',handle)
flops = 0._dbl
nx1 = SIZE ( nhc % nvt ( :, 1 ) )
nx2 = SIZE ( nhc % nvt ( 1, : ) )
! get force on first thermostat for all the chains in the system.
nhc % v_scale = 1.0_dbl
nhc % s_kin = 0.0_dbl
ncoef = SIZE ( ao % cr)
DO i = 1, ncoef
nhc % p_kin(1,i) % point = nhc % p_kin(1,i) % point + &
ao % mass ( i ) * ao % cr ( i ) * ao % cr ( i )
END DO
! force on the first bead in every thermostat chain (f1=sum(m*v**2)-ncoef*KbT)
DO n = 1, nhc % num_nhc
nhc % nvt ( 1, n ) % f = ( nhc % s_kin ( n ) - nhc % nvt ( 1, n ) % nkt )/ &
nhc % nvt ( 1, n ) % mass
END DO
! perform multiple time stepping using Yoshida
NCLOOP: DO inc = 1, nhc % nc
YOSH: DO iyosh = 1, nhc % nyosh
!
! update velocity on the last thermostat in the chain ! O1
nhc % nvt ( nhc % nhc_len, : ) % v = nhc % nvt (nhc % nhc_len, : ) % v + &
nhc % nvt (nhc % nhc_len, : ) % f * 0.25_dbl * nhc % dt_yosh ( iyosh )
!
! update velocity of other thermostats on chain (from nhc_len-1 to 1) ! O2
DO n = 1, nhc % num_nhc
DO inhc = nhc % nhc_len - 1, 1, -1
scale = exp (-0.125_dbl * nhc % nvt ( inhc+1, n )% v * &
nhc % dt_yosh ( iyosh ) )
nhc % nvt ( inhc, n ) % v = nhc % nvt ( inhc, n ) % v * scale ! scale
nhc % nvt ( inhc, n ) % v = nhc % nvt ( inhc, n ) % v + &
nhc % nvt ( inhc, n ) % f * 0.25_dbl * nhc % dt_yosh ( iyosh ) ! shift
nhc % nvt ( inhc, n ) % v = nhc % nvt ( inhc , n ) % v * scale ! scale
END DO
END DO
! the core of the operator ----- START------
! update nhc positions
nhc % nvt ( :, : ) % eta = nhc % nvt ( :, : ) % eta + &
0.5_dbl * nhc % nvt ( :, : ) % v * nhc % dt_yosh ( iyosh )
! now accumulate the scale factor for particle velocities
DO n = 1, nhc % num_nhc
nhc % v_scale ( n ) = nhc % v_scale ( n ) * &
exp ( -0.5_dbl * nhc % dt_yosh ( &
iyosh ) * nhc % nvt ( 1, n ) % v )
END DO
! the core of the operator ------ END ------
! update the force on first thermostat again (since coefficient velocities
! have changed)
DO n = 1, nhc % num_nhc
nhc % nvt ( 1, n ) % f = ( nhc % s_kin ( n ) * nhc % v_scale ( n ) &
* nhc % v_scale ( n ) - nhc % nvt ( 1, n ) % nkt ) &
/nhc % nvt ( 1, n ) % mass
END DO
! update velocity of other thermostats on chain (from 1 to nhc_len-1) ! O2
DO inhc = 1, nhc % nhc_len - 1
DO n = 1, nhc % num_nhc
scale = exp( -0.125_dbl *nhc % nvt ( inhc+1, n ) % v * &
nhc % dt_yosh ( iyosh ) )
nhc % nvt ( inhc, n ) % v = nhc % nvt( inhc, n ) % v * scale ! scale
nhc % nvt ( inhc, n ) % v = nhc % nvt( inhc, n ) % v + &
nhc % nvt ( inhc, n ) % f * 0.25_dbl * nhc % dt_yosh ( iyosh ) ! shift
nhc % nvt ( inhc, n ) % v = nhc % nvt ( inhc, n ) % v * scale ! scale
END DO
! updating the forces on all the thermostats
DO n = 1, nhc % num_nhc
nhc % nvt ( inhc+1, n ) % f = &
(nhc % nvt ( inhc, n ) % mass * nhc % nvt ( inhc, n ) % v &
*nhc % nvt ( inhc, n ) % v-nhc % nvt ( inhc + 1, n ) % nkt ) / &
nhc % nvt ( inhc + 1, n ) % mass
END DO
END DO
! update velocity on last thermostat ! O1
nhc % nvt ( nhc % nhc_len, : ) % v = &
nhc % nvt ( nhc % nhc_len, : ) % v + &
nhc % nvt ( nhc % nhc_len, : ) % f * &
0.25_dbl * nhc % dt_yosh ( iyosh )
END DO YOSH
END DO NCLOOP
! now scale the coefficients velocities
DO i = 1, ncoef
ao % cr ( i ) = ao % cr ( i ) * nhc % p_scale ( 1, i ) % point
END DO
!
CALL timestop(flops,handle)
END SUBROUTINE lnhc_ao
!******************************************************************************
!!****** nose/lnhc_barostat [1.0] *
!!
@ -380,10 +505,10 @@ SUBROUTINE lnhc_particles ( nhc, pnode, group )
END SUBROUTINE lnhc_particles
!******************************************************************************
!!****** nose/initial_nhc [1.0] *
!!****** nose/initialize_nhc [1.0] *
!!
!! NAME
!! initial_nhc
!! initialize_nhc
!!
!! FUNCTION
!! Initializes the NHC velocities to the Maxwellian distribution
@ -399,7 +524,7 @@ END SUBROUTINE lnhc_particles
!! SOURCE
!******************************************************************************
SUBROUTINE initial_nhc(nhc,tau_nhc,temp_ext, globenv )
SUBROUTINE initialize_nhc(nhc,tau_nhc,temp_ext, globenv )
IMPLICIT NONE
TYPE (lnhc_parameters_type ), INTENT ( INOUT ) :: nhc
TYPE ( global_environment_type ), INTENT ( INOUT ) :: globenv
@ -453,14 +578,14 @@ SUBROUTINE initial_nhc(nhc,tau_nhc,temp_ext, globenv )
END DO
END DO
END SUBROUTINE initial_nhc
END SUBROUTINE initialize_nhc
!!*****
!******************************************************************************
!!****** nose/initial_npt [1.0] *
!!****** nose/initialize_npt [1.0] *
!!
!! NAME
!! initial_npt
!! initialize_npt
!!
!! FUNCTION
!! Initializes the NHC velocities to the Maxwellian distribution
@ -476,7 +601,7 @@ END SUBROUTINE initial_nhc
!! SOURCE
!******************************************************************************
SUBROUTINE initial_npt ( npt, tau_cell, temp_ext, deth, nfree, ensemble, &
SUBROUTINE initialize_npt ( npt, tau_cell, temp_ext, deth, nfree, ensemble, &
globenv )
IMPLICIT NONE
@ -547,7 +672,7 @@ SUBROUTINE initial_npt ( npt, tau_cell, temp_ext, deth, nfree, ensemble, &
! initialize the position
IF ( ensemble == 'NPT_I' ) npt ( :, : ) % eps = LOG ( deth ) / 3.0_dbl
END SUBROUTINE initial_npt
END SUBROUTINE initialize_npt
!!*****
!******************************************************************************

4
src/pair_potential.F
View file

@ -145,7 +145,7 @@ FUNCTION ener_coul ( pot, qi, qj, r, ewald_param )
!------------------------------------------------------------------------------
SELECT CASE( ewald_param % ewald_type )
CASE( 'EWALD', 'PME' )
CASE( 'EWALD', 'PME', 'SPME' )
ener_coul = ifourpi / ewald_param % eps0 * qi * qj &
* erfc ( ewald_param % alpha * r ) / r - pot % e_cutoff_coul
CASE( 'NONE')
@ -274,7 +274,7 @@ SUBROUTINE potential_s ( rijsq, potparm, iatom, jatom, energy, fscalar )
energy = splineh ( spline_table ( spltab ( iatom, jatom ) ), rijsq, fscalar )
fscalar = -2.0_dbl * fscalar
END SUBROUTINE potential_s
!******************************************************************************

66
src/particle_lists.F
View file

@ -23,15 +23,17 @@ MODULE particle_lists
USE kinds, ONLY : dbl
USE molecule_types, ONLY : particle_node_type
USE particle_types, ONLY : particle_type
USE termination, ONLY : stop_memory
IMPLICIT NONE
PRIVATE
PUBLIC :: particle_list_type, particle_list_select
PUBLIC :: particle_list_type, particle_list_ewald, particle_list_grid
TYPE particle_list_type
INTEGER :: particle_index
INTEGER :: grid_id
REAL ( dbl ), DIMENSION ( : ), POINTER :: r
REAL ( dbl ), POINTER :: charge
TYPE ( particle_list_type ), POINTER :: next
@ -44,7 +46,7 @@ CONTAINS
!******************************************************************************
SUBROUTINE particle_list_select ( pnode, plist, nparts)
SUBROUTINE particle_list_ewald ( pnode, plist, nparts)
IMPLICIT NONE
! Arguments
@ -93,7 +95,65 @@ SUBROUTINE particle_list_select ( pnode, plist, nparts)
END IF
END DO
END SUBROUTINE particle_list_select
END SUBROUTINE particle_list_ewald
!****************************************************************************
SUBROUTINE particle_list_grid ( part, plist, atom_list, nparts, grid_tag )
IMPLICIT NONE
! Arguments
TYPE ( particle_type ), DIMENSION ( : ), INTENT ( INOUT ), TARGET :: part
TYPE ( particle_list_type ), POINTER :: plist
INTEGER, INTENT (IN), DIMENSION (:) :: ATOM_LIST
INTEGER , INTENT ( OUT ) :: nparts
INTEGER , INTENT ( IN ) :: grid_tag
! Locals
INTEGER :: npart, isos, part_index
TYPE ( particle_list_type ), POINTER :: plist_loc
REAL (dbl), TARGET :: charge
!------------------------------------------------------------------------------
! 'plist' will always point to the beginning of the list, 'plist_loc'
! to the current position
NULLIFY ( plist )
NULLIFY ( plist_loc )
nparts = 0
charge = 1.0_dbl
DO npart = 1, SIZE ( atom_list )
IF ( .NOT. ASSOCIATED ( plist_loc ) ) THEN
ALLOCATE ( plist, STAT = isos )
IF ( isos /= 0 ) &
CALL stop_memory ( "build_particle_list", "plist", 0 )
plist_loc => plist
NULLIFY ( plist_loc % next )
ELSE
ALLOCATE ( plist_loc % next, STAT = isos )
IF ( isos /= 0 ) &
CALL stop_memory ( "build_particle_list", "next", 0 )
plist_loc => plist_loc % next
NULLIFY ( plist_loc % next )
END IF
nparts = nparts + 1
part_index = atom_list (npart)
part ( part_index ) % grid_id = grid_tag
plist_loc % particle_index = atom_list (npart)
plist_loc % grid_id = grid_tag
plist_loc % r => part ( part_index ) % r
plist_loc % charge => charge
END DO
END SUBROUTINE particle_list_grid
!******************************************************************************

4
src/particle_types.F
View file

@ -36,11 +36,15 @@ MODULE particle_types
TYPE particle_type
INTEGER :: iatom
INTEGER :: grid_id
INTEGER :: node
REAL ( dbl ), DIMENSION (3) :: r
REAL ( dbl ), DIMENSION (3) :: f
REAL ( dbl ), DIMENSION (3) :: v
TYPE (particle_prop_type ), POINTER :: prop
TYPE ( cell_neighborlist_type ), POINTER :: cell
TYPE ( cell_neighborlist_type ), POINTER :: cell_ol
INTEGER, DIMENSION(:), POINTER :: coef_list
END TYPE particle_type
END MODULE particle_types

12
src/pme.F
View file

@ -1,4 +1,6 @@
!-----------------------------------------------------------------------------!
!-----------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
@ -25,13 +27,12 @@ MODULE pme
PW_COMPLEXDATA3D, PW_REALDATA3D, &
PW_REALSPACE, PW_COMPLEXDATA1D, PW_RECIPROCALSPACE
USE dgs, ONLY : dg_get_strucfac, dg_get_patch, dg_sum_patch, &
dg_sum_patch_force
dg_sum_patch_force_3d
USE dg_types, ONLY : dg_type
USE ewald_parameters_types, ONLY : ewald_parameters_type
USE hartree, ONLY : calculate_hartree
USE kinds, ONLY: dbl, sgl, dp_size, sp_size
USE mathconstants, ONLY : fourpi, pi
USE particle_lists, ONLY : particle_list_type
USE particle_types, ONLY : particle_type
USE pw_grid_types, ONLY : pw_grid_type
USE pws, ONLY : pw_poisson_solver, init_pw_poisson_solver, pw_integral_a2b
@ -164,7 +165,6 @@ SUBROUTINE pme_evaluate ( dg, part, fg_coulomb, vg_coulomb, pv_g, box, &
CALL calculate_hartree ( rhob_r, vg_coulomb, phi_r, dphi_g, h_stress )
CALL coeff_deallocate ( phi_r )
vg_coulomb = vg_coulomb * inv_fourpi_eps0
!---------- END OF ELECTROSTATIC CALCULATION --------
@ -184,6 +184,8 @@ SUBROUTINE pme_evaluate ( dg, part, fg_coulomb, vg_coulomb, pv_g, box, &
!--------END OF STRESS TENSOR CALCULATION -----------
CALL coeff_deallocate ( phi_r )
CALL rs_grid_setup ( drpot, grid_b , grid_s % npts)
CALL rs_grid_allocate ( drpot ( 1 : 3 ) )
@ -212,12 +214,12 @@ SUBROUTINE pme_evaluate ( dg, part, fg_coulomb, vg_coulomb, pv_g, box, &
CALL get_patch_again ( dg, part, exp_igr, p1, p2, rhos1, rhos2 )
! sum boxes on real space grids (big box)
CALL dg_sum_patch_force ( drpot, rhos1, exp_igr % centre ( :, p1 ), fat )
CALL dg_sum_patch_force_3d ( drpot, rhos1, exp_igr % centre ( :, p1 ), fat )
fg_coulomb ( 1, p1 ) = fg_coulomb ( 1, p1 ) - fat ( 1 ) * dvols
fg_coulomb ( 2, p1 ) = fg_coulomb ( 2, p1 ) - fat ( 2 ) * dvols
fg_coulomb ( 3, p1 ) = fg_coulomb ( 3, p1 ) - fat ( 3 ) * dvols
IF ( p2 /= 0 ) THEN
CALL dg_sum_patch_force ( drpot, rhos2, exp_igr % centre ( :, p2 ), fat )
CALL dg_sum_patch_force_3d ( drpot, rhos2, exp_igr % centre ( :, p2 ), fat )
fg_coulomb ( 1, p2 ) = fg_coulomb ( 1, p2 ) - fat ( 1 ) * dvols
fg_coulomb ( 2, p2 ) = fg_coulomb ( 2, p2 ) - fat ( 2 ) * dvols
fg_coulomb ( 3, p2 ) = fg_coulomb ( 3, p2 ) - fat ( 3 ) * dvols

183
src/pol_coefs.F Normal file
View file

@ -0,0 +1,183 @@
!-----------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!***** cp2k/pol_coefs [1.0] *
!!
!! NAME
!! pol_coefs
!!
!! FUNCTION
!!
!! AUTHOR
!! gloria
!!
!! MODIFICATION HISTORY
!!
!! NOTES
!!
!!
!! SOURCE
!******************************************************************************
MODULE pol_coefs
USE atomic_kinds, ONLY : kind_info_type
USE ao_types, ONLY : ao_type, ao_allocate, ao_zero
USE coefficient_types, ONLY : coeff_type
USE kinds, ONLY : dbl
USE particle_types, ONLY : particle_type
USE molecule_types , ONLY : molecule_type, read_basis_type
USE termination, ONLY : stop_memory
IMPLICIT NONE
PUBLIC :: allocate_pol_coefs, initialize_pol_coefs
PRIVATE
CONTAINS
!******************************************************************************
SUBROUTINE allocate_pol_coefs(ncoef, c_pos, c_vel, c_force)
IMPLICIT NONE
INTEGER, INTENT (in) :: ncoef
TYPE ( coeff_type ), INTENT(OUT) :: c_pos
TYPE ( coeff_type ), INTENT(OUT) :: c_vel
TYPE ( coeff_type ), INTENT(OUT) :: c_force
! locals
INTEGER, PARAMETER :: ATOMICORBITALS = 402
LOGICAL :: mass_flag, grid_flag
! allocate cpos
c_pos % in_use = ATOMICORBITALS
mass_flag = .TRUE.
grid_flag = .TRUE.
CALL ao_allocate ( c_pos % ao, ncoef, mass_flag, grid_flag )
CALL ao_zero ( c_pos % ao)
! allocate cvel
c_vel % in_use = ATOMICORBITALS
mass_flag = .FALSE.
grid_flag = .FALSE.
CALL ao_allocate ( c_vel % ao, ncoef, mass_flag, grid_flag )
CALL ao_zero ( c_vel % ao)
! allocate cforce
c_force % in_use = ATOMICORBITALS
mass_flag = .FALSE.
grid_flag = .FALSE.
CALL ao_allocate ( c_force % ao, ncoef, mass_flag, grid_flag )
CALL ao_zero ( c_force % ao)
END SUBROUTINE allocate_pol_coefs
!***************************************************************************
SUBROUTINE initialize_pol_coefs ( rtype , ki, part, c_pos, c_vel, c_force)
IMPLICIT NONE
CHARACTER (LEN = 12), INTENT (in) :: rtype
TYPE (kind_info_type), POINTER, DIMENSION(:) :: ki
TYPE (particle_type), intent(in), DIMENSION(:) :: part
TYPE ( coeff_type ), INTENT(OUT) :: c_pos
TYPE ( coeff_type ), INTENT(OUT) :: c_vel
TYPE ( coeff_type ), INTENT(OUT) :: c_force
! locals
INTEGER :: ncoef, icoef
LOGICAL :: kind_flag
ncoef = size (c_pos%ao%cr(:))
SELECT CASE (rtype)
CASE ("INIT")
!! CALL initialize_coeff_velocities (c_vel, stuff)!!
c_pos % ao % cr ( : ) = 1.0_dbl
CASE ("POS")
!! CALL initialize_coeff_velocities (c_vel, stuff)!!
OPEN (unit=666,file='coefficent.rst',status='unknown')
DO icoef = 1, ncoef
READ (666,*) c_pos % ao % cr (icoef)
END DO
CLOSE (666)
CASE ("POSVEL")
OPEN (unit=666,file='coefficent.rst',status='unknown')
DO icoef = 1, ncoef
READ (666,*) c_pos % ao % cr (icoef), c_vel % ao % cr (icoef)
END DO
CLOSE (666)
END SELECT
! fills the remaining information...
kind_flag = .TRUE.
CALL get_kind_and_part_index ( c_pos % ao, ki, part, kind_flag)
kind_flag = .FALSE.
CALL get_kind_and_part_index ( c_vel % ao, ki, part, kind_flag)
kind_flag = .FALSE.
CALL get_kind_and_part_index ( c_force % ao, ki, part, kind_flag)
END SUBROUTINE initialize_pol_coefs
!------------------------------------------------------------------------------!
SUBROUTINE get_kind_and_part_index (ao, ki, part, flag)
IMPLICIT NONE
TYPE (ao_type), INTENT (INOUT) :: ao
TYPE (kind_info_type), POINTER, DIMENSION(:) :: ki
TYPE (particle_type), intent(in), DIMENSION(:) :: part
LOGICAL, intent (in) :: flag
! locals
INTEGER :: ncoef, ikind, ipart, ii, icoef, icgf, ncgf, ios
TYPE (kind_info_type), POINTER :: kinfo
ncoef = size (ao % cr)
DO ikind = 1 , size (ki)
DO ii = 1, ki(ikind)%natom
ipart = ki(ikind)%atom_list(ii)
ncgf = ki(ikind) % orb_basis_set % ncgf
DO icgf = 1, ncgf
icoef = part ( ipart ) % coef_list( icgf )
ao % coef_to_basis_set ( icoef ) = ikind
ao % coef_to_part ( icoef ) = ipart
IF (flag) THEN
ALLOCATE (ao % kind_info (ncoef), STAT = ios)
kinfo => ao % kind_info ( icoef )
kinfo => ki (ikind)
END IF
ENDDO
END DO
END DO
END SUBROUTINE get_kind_and_part_index
!------------------------------------------------------------------------------!
END MODULE pol_coefs
!------------------------------------------------------------------------------!

1121
src/pol_electrostatics.F Normal file
View file

File diff suppressed because it is too large Load diff

254
src/pol_force.F Normal file
View file

@ -0,0 +1,254 @@
!-----------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C ) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/fist_force [1.0] *
!!
!! NAME
!! pol_force
!!
!! FUNCTION
!!
!! AUTHOR
!! gloria
!!
!! SOURCE
!******************************************************************************
MODULE pol_force
USE atomic_kinds, ONLY : kind_info_type
USE dg_types, ONLY : dg_type
USE coefficient_types, ONLY : coeff_type
USE ewald_parameters_types, ONLY : ewald_parameters_type
USE fist_intra_force, ONLY : force_intra_control
USE fist_nonbond_force, ONLY : force_nonbond, bonded_correct_gaussian
USE global_types, ONLY : global_environment_type
USE greens_fn, ONLY : pw_green_fn_init
USE kinds, ONLY : dbl
USE linklist_control, ONLY : list_control
USE mathconstants, ONLY : pi, zero
USE md, ONLY : thermodynamic_type
USE molecule_types, ONLY : molecule_structure_type, particle_node_type, &
linklist_neighbor, linklist_exclusion, linklist_bonds, &
linklist_bends, linklist_torsions
USE message_passing, ONLY : mp_sum
USE pair_potential, ONLY : potentialparm_type
USE pol_electrostatics, ONLY : electrostatics
USE particle_types, ONLY : particle_type
USE pol_grids, ONLY : pol_grids_allocate, pol_grids_count, &
pol_grids_map, pol_grids_initialize, assign_coefs_to_grids
USE pme, ONLY: pme_evaluate
USE pw_grid_types, ONLY : pw_grid_type
USE pw_grids, ONLY : pw_grid_change
USE simulation_cell, ONLY : cell_type, get_hinv
USE spme, ONLY : spme_evaluate
USE termination, ONLY : stop_program, stop_memory
USE timings, ONLY : timeset, timestop
USE linklist_types, ONLY : linklist_internal_data_type
PRIVATE
PUBLIC :: force_control, debug_variables_type
TYPE debug_variables_type
REAL ( dbl ) :: pot_nonbond, pot_g, pot_bond, pot_bend, pot_torsion
REAL ( dbl ), DIMENSION ( :, : ), POINTER :: &
f_nonbond, f_g, f_bond, f_bend, f_torsion
REAL ( dbl ), DIMENSION ( 3, 3 ) :: pv_nonbond, pv_g, pv_bond, &
pv_bend, pv_torsion
END TYPE debug_variables_type
!!*****
!******************************************************************************
CONTAINS
!******************************************************************************
!!****** pol_force/force_control [1.0] *
!!
!! NAME
!! force_control
!!
!! SOURCE
!******************************************************************************
SUBROUTINE force_control ( molecule, pnode, part, box, box_ref, &
drho_basis_info, rho0_basis_info, coef_pos, coef_vel, coef_force, &
thermo, potparm, ewald_param, ensemble, fc_global, lldata, debug)
IMPLICIT NONE
! Arguments
TYPE ( particle_type ), DIMENSION ( : ), INTENT ( INOUT ) :: part
TYPE ( particle_node_type ), DIMENSION ( : ), INTENT ( INOUT ) :: pnode
TYPE ( molecule_structure_type ), DIMENSION ( : ), INTENT ( IN ) :: molecule
TYPE ( cell_type ), INTENT ( INOUT ) :: box
TYPE ( cell_type ), INTENT ( IN ) :: box_ref
TYPE ( kind_info_type ), POINTER, DIMENSION ( : ) :: drho_basis_info
TYPE ( kind_info_type ), POINTER, DIMENSION ( : ) :: rho0_basis_info
TYPE ( coeff_type ), INTENT (inout) :: coef_pos
TYPE ( coeff_type ), INTENT (inout) :: coef_vel
TYPE ( coeff_type ), INTENT (inout) :: coef_force
TYPE ( thermodynamic_type ), INTENT ( INOUT ) :: thermo
TYPE ( potentialparm_type ), DIMENSION ( :,: ), INTENT ( IN ) :: potparm
TYPE ( ewald_parameters_type ), INTENT ( INOUT ) :: ewald_param
CHARACTER ( LEN = * ), INTENT ( IN ) :: ensemble
TYPE ( global_environment_type ), INTENT ( IN ) :: fc_global
TYPE ( linklist_internal_data_type), DIMENSION (:), INTENT (INOUT) :: lldata
TYPE ( debug_variables_type ), INTENT ( OUT ), OPTIONAL :: debug
! Locals
INTEGER :: id, i, ii, natoms, nnodes, handle, isos, ncoefs
INTEGER :: ng, igrid, ng_coef, ng_atom, nt_coef, nt_atom, iref
REAL ( dbl ) :: pot_nonbond
REAL ( dbl ), DIMENSION ( :,: ), ALLOCATABLE, SAVE :: f_nonbond
REAL ( dbl ), DIMENSION ( :,: ), ALLOCATABLE :: fc_part
REAL ( dbl ), DIMENSION ( :,: ), ALLOCATABLE :: fo_part
REAL ( dbl ), DIMENSION (:), ALLOCATABLE :: fc_coef
REAL ( dbl ), DIMENSION (:), ALLOCATABLE :: fo_coef
REAL ( dbl ), DIMENSION ( :,: ), ALLOCATABLE :: f_total
REAL ( dbl ), DIMENSION ( 3,3 ) :: pv_nonbond
TYPE ( pw_grid_type ),POINTER, DIMENSION (:), SAVE :: big_grid, small_grid
TYPE ( dg_type ), POINTER, DIMENSION (:), SAVE :: dg_coef
TYPE ( dg_type ), POINTER, DIMENSION (:), SAVE :: dg_atom
LOGICAL :: first_time
!------------------------------------------------------------------------------
CALL timeset ( 'FORCE','I',' ',handle )
nnodes = SIZE ( pnode )
natoms = SIZE ( part )
ncoefs = size ( coef_pos % ao % cr )
isos = 0
first_time = .NOT. ALLOCATED ( f_nonbond )
IF ( .NOT. ALLOCATED ( f_nonbond ) ) THEN
ALLOCATE ( f_nonbond ( 3,natoms ), STAT = isos )
IF ( isos /= 0 ) &
CALL stop_memory ( 'force_control', 'f_nonbond', 3 * natoms )
ELSE IF ( SIZE ( f_nonbond ( 1, : ) ) < natoms ) THEN
DEALLOCATE ( f_nonbond, STAT = isos )
IF ( isos /= 0 ) CALL stop_memory ( 'force_control', 'f_nonbond' )
ALLOCATE ( f_nonbond ( 3,natoms ), STAT = isos )
IF ( isos /= 0 ) &
CALL stop_memory ( 'force_control', 'f_nonbond', 3 * natoms )
END IF
IF ( first_time ) THEN
CALL pol_grids_count ( drho_basis_info, rho0_basis_info, ng_coef, ng_atom, &
nt_coef, nt_atom )
ng = ng_coef + ng_atom
CALL pol_grids_allocate (dg_coef, dg_atom, small_grid, big_grid, &
nt_coef, nt_atom, ng)
igrid = 0
CALL pol_grids_map (dg_coef, drho_basis_info, igrid)
CALL pol_grids_map (dg_atom, rho0_basis_info, igrid)
CALL assign_coefs_to_grids (coef_pos % ao, part, drho_basis_info, nt_coef)
CALL pol_grids_initialize ( dg_coef, dg_atom, part, pnode, coef_pos % ao, &
drho_basis_info, rho0_basis_info, fc_global, box, thermo, &
small_grid, big_grid, iref )
CALL pw_green_fn_init ( box % green, big_grid ( iref ) )
END IF
!
! first check with list_control to update neighbor lists
!
CALL list_control ( lldata (1), pnode, part, box )
!
! check to update the overlap neighbor list
!
CALL list_control ( lldata (2), pnode, part, box, drho_basis_info )
!
! initial force, energy and pressure tensor arrays
!
DO i = 1, natoms
part ( i ) % f ( 1 ) = 0.0_dbl
part ( i ) % f ( 2 ) = 0.0_dbl
part ( i ) % f ( 3 ) = 0.0_dbl
END DO
thermo % pot = 0.0_dbl
thermo % pv = 0.0_dbl
IF ( .NOT. ALLOCATED ( fc_part ) ) &
ALLOCATE ( fc_part ( 3,nnodes ), STAT=isos )
IF ( isos /= 0 ) &
CALL stop_memory ( 'force_control', 'fc_part', 3 * nnodes )
IF ( .NOT. ALLOCATED ( fc_coef ) ) &
ALLOCATE ( fc_coef ( ncoefs ), STAT=isos )
IF ( isos /= 0 ) &
CALL stop_memory ( 'force_control', 'fc_coef', ncoefs )
! get electrostatics forces
CALL electrostatics (dg_atom, dg_coef, thermo, part, coef_pos, box, &
small_grid, big_grid, ewald_param, iref, fc_coef, fc_part )
!
! get real-space non-bonded forces:
!
f_nonbond = zero
CALL force_nonbond ( ewald_param,pnode,box,potparm, &
pot_nonbond,f_nonbond,pv_nonbond )
!! CALL pol_overlap_force (STUFF)
!
! add up all the potential energies
!
! add up all the forces
! nonbonded forces might be calculated for atoms not on this node
! ewald forces are strictly local -> sum only over pnode
! We first sum the forces in f_nonbond, this allows for a more efficient
! global sum in the parallel code and in the end copy them back to part
isos = 0
IF ( .NOT. ALLOCATED ( f_total ) ) ALLOCATE ( f_total ( 3,natoms ), STAT=isos )
IF ( isos /= 0 ) CALL stop_memory ( 'force_control', 'f_total', 3 * natoms )
DO i = 1, natoms
f_total ( 1, i ) = part ( i ) % f ( 1 ) + f_nonbond ( 1, i )
f_total ( 2, i ) = part ( i ) % f ( 2 ) + f_nonbond ( 2, i )
f_total ( 3, i ) = part ( i ) % f ( 3 ) + f_nonbond ( 3, i )
END DO
DO i = 1, nnodes
ii = pnode ( i ) % p % iatom
f_total ( 1, ii ) = f_total ( 1, ii ) + fc_part ( 1, i )
f_total ( 2, ii ) = f_total ( 2, ii ) + fc_part ( 2, i )
f_total ( 3, ii ) = f_total ( 3, ii ) + fc_part ( 3, i )
END DO
! deallocating all local variables
isos = 0
IF ( ALLOCATED ( fc_part ) ) DEALLOCATE ( fc_part, STAT = isos )
IF ( isos /= 0 ) CALL stop_memory ( 'force_control', 'fc_part' )
isos = 0
IF ( ALLOCATED ( f_total ) ) DEALLOCATE ( f_total, STAT = isos )
IF ( isos /= 0 ) CALL stop_memory ( 'force_control', 'f_total' )
CALL timestop ( zero, handle )
END SUBROUTINE force_control
!!*****
!******************************************************************************
END MODULE pol_force
!******************************************************************************

576
src/pol_grids.F Normal file
View file

@ -0,0 +1,576 @@
!-----------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/pol_grids [1.0] *
!!
!! NAME
!! pol_grids
!!
!! FUNCTION
!!
!! AUTHOR
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE pol_grids
USE atomic_kinds, ONLY : kind_info_type
USE ao_types, ONLY : ao_type
USE basis_set_types, ONLY : gto_basis_set_type
USE global_types, ONLY : global_environment_type
USE dg_types, ONLY : dg_type
USE dgs, ONLY : dg_pol_grid_setup, dg_find_cutoff
USE dg_rho0s, ONLY : dg_rho0_type, dg_rho0_setup
USE kinds, ONLY : dbl
USE mathconstants, ONLY : pi, zero
USE md, ONLY : thermodynamic_type
USE molecule_types, ONLY : particle_node_type
USE message_passing, ONLY : mp_sum
USE particle_types, ONLY : particle_type
USE particle_lists, ONLY : particle_list_type, particle_list_grid
USE coefficient_lists, ONLY : coef_list_type, coef_list_grid
USE pw_grid_types, ONLY : pw_grid_type, HALFSPACE
USE simulation_cell, ONLY : cell_type
USE structure_factor_types, ONLY : structure_factor_type
USE structure_factors, ONLY : structure_factor_evaluate, &
structure_factor_allocate, &
structure_factor_deallocate
USE termination, ONLY : stop_memory, stop_program
USE timings, ONLY : timeset, timestop
USE util, ONLY : matvec_3x3
IMPLICIT NONE
PRIVATE
PUBLIC :: pol_grids_allocate, pol_grids_count, pol_grids_initialize
PUBLIC :: pol_grids_map, assign_coefs_to_grids
!!*****
!******************************************************************************
CONTAINS
!******************************************************************************
!*****
!******************************************************************************
!!****** pol_grids/pol_grids_allocate [1.0] *
!!
!! NAME
!! pol_grids_allocate
!!
!! FUNCTION
!!
!! AUTHOR
!! gloria
!!
!! SOURCE
!******************************************************************************
SUBROUTINE pol_grids_allocate ( dg_coef, dg_atom, small_grid, big_grid, &
ntc, nta, ng )
IMPLICIT NONE
! Arguments
TYPE ( dg_type ), POINTER, DIMENSION (:) :: dg_coef
TYPE ( dg_type ), POINTER, DIMENSION (:) :: dg_atom
TYPE ( pw_grid_type ), POINTER, DIMENSION (:) :: small_grid
TYPE ( pw_grid_type ), POINTER, DIMENSION (:) :: big_grid
INTEGER, intent(IN) :: ntc
INTEGER, intent(IN) :: nta
INTEGER, intent(IN) :: ng
! Locals
INTEGER :: ios
!---------------------------------------------
NULLIFY (dg_coef, dg_atom)
NULLIFY (small_grid, big_grid)
ALLOCATE (dg_coef(ntc), STAT = ios)
IF (ios /= 0 ) THEN
CALL stop_memory ('allocate_pol_grids', 'dg_coef', ntc)
END IF
ALLOCATE (dg_atom(nta), STAT = ios)
IF (ios /= 0 ) THEN
CALL stop_memory ('allocate_pol_grids', 'dg_atom', nta)
END IF
ALLOCATE (small_grid(ng), STAT = ios)
IF (ios /= 0 ) THEN
CALL stop_memory ('allocate_pol_grids', 'small_grid', ng)
END IF
ALLOCATE (big_grid(ng), STAT = ios)
IF (ios /= 0 ) THEN
CALL stop_memory ('allocate_pol_grids', 'big_grid', ng)
END IF
END SUBROUTINE pol_grids_allocate
!******************************************************************************
!*****
!******************************************************************************
!!****** pol_grids/pol_grids_count [1.0] *
!!
!! NAME
!! pol_grids_count
!!
!! FUNCTION
!!
!! AUTHOR
!! gloria
!!
!! SOURCE
!******************************************************************************
SUBROUTINE pol_grids_count (drho_basis_info, rho0_basis_info, ngc, nga, ntc, nta)
IMPLICIT NONE
TYPE ( kind_info_type ), INTENT(IN), DIMENSION ( : ) :: drho_basis_info
TYPE ( kind_info_type ), INTENT(IN), DIMENSION ( : ) :: rho0_basis_info
INTEGER, INTENT(out) :: ngc
INTEGER, INTENT(out) :: nga
INTEGER, INTENT(out) :: ntc
INTEGER, INTENT(out) :: nta
! locals
INTEGER :: i, iset, nshell_set, ishell, nkinda, nkindc
nkindc = size(drho_basis_info)
nkinda = size(rho0_basis_info)
nshell_set = 0
nga = 0
ngc = 0
nta = 0
ntc = 0
! coeffs loop
DO i=1, nkindc
DO iset = 1, drho_basis_info(i) % orb_basis_set% nset
nshell_set = drho_basis_info(i)%orb_basis_set % nshell(iset)
ngc = ngc + nshell_set
DO ishell = 1, nshell_set
ntc = ntc +(2*drho_basis_info(i)% orb_basis_set % l(ishell,iset)+1)
ENDDO
ENDDO
ENDDO
nshell_set = 0
! atom loop
DO i=1, nkinda
DO iset = 1, rho0_basis_info(i) % orb_basis_set%nset
nshell_set = rho0_basis_info(i)%orb_basis_set % nshell(iset)
nga = nga + nshell_set
DO ishell = 1, nshell_set
nta = nta +(2*rho0_basis_info(i)% orb_basis_set % l(ishell,iset)+1)
ENDDO
ENDDO
ENDDO
END SUBROUTINE pol_grids_count
!******************************************************************************
!*****
!******************************************************************************
!!****** pol_grids/pol_grids_map [1.0] *
!!
!! NAME
!! pol_grids_map
!!
!! FUNCTION
!!
!! AUTHOR
!! gloria
!!
!! SOURCE
!******************************************************************************
SUBROUTINE pol_grids_map ( dg, ki, igrid)
IMPLICIT NONE
TYPE ( dg_type ), POINTER, DIMENSION (:) :: dg
TYPE ( kind_info_type ), POINTER, DIMENSION ( : ) :: ki
INTEGER, INTENT(INOUT) :: igrid
! locals
INTEGER :: igt, ikind, iset, ishell, nshell, ncc, icc, ios, ip
TYPE ( gto_basis_set_type ), POINTER :: gto_info
igt = 0
DO ikind=1,SIZE(ki)
gto_info => ki(ikind)%orb_basis_set
DO iset = 1, gto_info%nset
nshell = gto_info%nshell(iset)
DO ishell=1, nshell
igrid = igrid+1
SELECT CASE(gto_info%l(ishell,iset))
CASE(0)
igt=igt+1
dg(igt) % grid_index = igrid
dg(igt)% dg_rho0 % type = "S"
dg(igt)% dg_rho0 % grid = igrid
dg(igt)% dg_rho0 % kind = ikind
dg(igt)% dg_rho0 % cutoff_radius = gto_info%set_radius(iset)
dg(igt)% nparts = size (ki(ikind)%atom_list)
ncc=SIZE(gto_info % gcc(:,ishell,iset))
ALLOCATE(dg(igt)%dg_rho0%gcc(ncc), STAT=ios)
IF (ios /= 0 ) THEN
CALL stop_memory ('pol_grids_map', 'gcc', ncc)
END IF
ALLOCATE(dg(igt)%dg_rho0%zet(ncc), STAT=ios)
IF (ios /= 0 ) THEN
CALL stop_memory ('pol_grids_map', 'zet', ncc)
END IF
DO icc=1, ncc
dg(igt)% dg_rho0 % zet(icc) = gto_info%zet(icc,iset)
dg(igt)% dg_rho0 % gcc(icc) = gto_info%gcc(icc,ishell,iset)
ENDDO
CASE(1)
DO ip=1,3
igt=igt+1
SELECT CASE(ip)
CASE(1)
dg(igt)% dg_rho0% type ="PX"
CASE(2)
dg(igt)% dg_rho0% type ="PY"
CASE(3)
dg(igt)% dg_rho0% type ="PZ"
END SELECT
dg(igt) % grid_index = igrid
dg(igt)% dg_rho0% grid = igrid
dg(igt)% dg_rho0% kind = ikind
dg(igt)% dg_rho0 % cutoff_radius = gto_info%set_radius(iset)
dg(igt)% nparts = size (ki(ikind)%atom_list)
ncc=SIZE(gto_info % gcc(:,ishell,iset))
ALLOCATE(dg(igt)%dg_rho0%gcc(ncc), STAT=ios)
IF (ios /= 0 ) THEN
CALL stop_memory ('pol_grids_map', 'gcc', ncc)
END IF
ALLOCATE(dg(igt)%dg_rho0%zet(ncc), STAT=ios)
IF (ios /= 0 ) THEN
CALL stop_memory ('pol_grids_map', 'zet', ncc)
END IF
DO icc=1, ncc
dg(igt)% dg_rho0 % zet(icc) = gto_info%zet(icc,iset)
dg(igt)% dg_rho0 % gcc(icc) = gto_info%gcc(icc,ishell,iset)
ENDDO
END DO
END SELECT
ENDDO
ENDDO
ENDDO
END SUBROUTINE pol_grids_map
!******************************************************************************
!*****
!!****** pol_grids/pol_grids_initialize [1.0] *
!!
!! NAME
!! pol_grids_initialize
!!
!! FUNCTION
!!
!! AUTHOR
!! gloria
!!
!! SOURCE
!******************************************************************************
SUBROUTINE pol_grids_initialize ( dg_coef, dg_atom, part, pnode, ao, &
drho_basis_info, rho0_basis_info, pnode_grp, box, &
thermo, small_grid, big_grid, iref )
IMPLICIT NONE
! Arguments
TYPE ( dg_type ), INTENT ( OUT ), DIMENSION (:) :: dg_coef
TYPE ( dg_type ), INTENT ( OUT ), DIMENSION (:) :: dg_atom
TYPE ( particle_type ), DIMENSION ( : ), INTENT ( INOUT ) :: part
TYPE ( particle_node_type ), DIMENSION ( : ), INTENT ( IN ) :: pnode
TYPE ( ao_type ), INTENT (inout) :: ao
TYPE ( kind_info_type ), POINTER, DIMENSION ( : ) :: drho_basis_info
TYPE ( kind_info_type ), POINTER, DIMENSION ( : ) :: rho0_basis_info
TYPE ( global_environment_type ), INTENT ( IN ) :: pnode_grp
TYPE ( cell_type ), INTENT ( IN ) :: box
TYPE ( thermodynamic_type ), INTENT ( INOUT ) :: thermo
TYPE ( pw_grid_type ), INTENT ( INOUT ), DIMENSION (:) :: small_grid
TYPE ( pw_grid_type ), INTENT ( INOUT ), DIMENSION (:) :: big_grid
INTEGER, INTENT (OUT) :: iref
! Locals
INTEGER :: natoms, npts_s ( 3 ), npts_max ( 3 ), iw
INTEGER :: igt, ikind, igrid, ig, igrid_old
REAL ( dbl ) :: qi, qj, i, cutoff, cutoff_radius, cutoff_max, &
max_cut0, max_cut
TYPE (kind_info_type), POINTER :: ki
!------------------------------------------------------------------------------
! fire up the reciprocal space and compute self interaction and
! term from the neutralizing background.
iw = pnode_grp % scr
! initialize igrid_old to negative number as to not
! be activated by accident
igrid_old = -1
! set up grids
! First: get the grid with max cutoff
iref = 0
max_cut0 = 1E-30_dbl
DO igt = 1, SIZE ( dg_coef )
ikind = dg_coef ( igt ) % dg_rho0 % kind
ki => drho_basis_info( ikind )
npts_s( : ) = ki % number_of_grid_points
igrid = dg_coef ( igt ) % dg_rho0 % grid
cutoff_radius = dg_coef ( igt ) % dg_rho0 % cutoff_radius
IF ( igrid == igrid_old ) CYCLE
CALL dg_find_cutoff ( box, npts_s, cutoff_radius, small_grid ( igrid ), &
big_grid ( igrid ), cutoff )
max_cut = MAX ( max_cut0, cutoff )
IF (max_cut /= max_cut0) THEN
iref = igrid
max_cut0 = max_cut
npts_max = npts_s
cutoff_max = cutoff_radius
END IF
igrid_old = igrid
ENDDO
DO igt = 1, SIZE ( dg_atom )
ikind = dg_atom ( igt ) % dg_rho0 % kind
ki => rho0_basis_info( ikind )
npts_s( : ) = ki % number_of_grid_points
igrid = dg_atom ( igt ) % dg_rho0 % grid
cutoff_radius = dg_atom ( igt ) % dg_rho0 % cutoff_radius
IF ( igrid == igrid_old ) CYCLE
CALL dg_find_cutoff ( box, npts_s, cutoff_radius, small_grid ( igrid ), &
big_grid ( igrid ), cutoff )
max_cut = MAX ( max_cut0, cutoff )
IF (max_cut /= max_cut0) THEN
iref = igrid
max_cut0 = max_cut
max_cut0 = max_cut
npts_max = npts_s
cutoff_max = cutoff_radius
END IF
igrid_old = igrid
ENDDO
! Second: setup the reference grid
CALL dg_pol_grid_setup ( box, npts_max, cutoff_max , &
small_grid ( iref ), big_grid ( iref ), pnode_grp )
! Third: setup rest of the grids w.r.t. the reference grid
! coef loop
DO igt = 1, SIZE ( dg_coef )
ikind = dg_coef ( igt ) % dg_rho0 % kind
ki => drho_basis_info( ikind )
npts_s( : ) = ki % number_of_grid_points
igrid = dg_coef ( igt ) % dg_rho0 % grid
cutoff_radius = dg_coef ( igt ) % dg_rho0 % cutoff_radius
CALL coef_list_grid ( ao, dg_coef ( igt ) % clist, part, &
dg_coef ( igt ) % ncoefs, igt, igrid )
IF ( igrid /= igrid_old ) THEN
IF ( igrid /= iref ) &
CALL dg_pol_grid_setup ( box, npts_s, cutoff_radius , &
small_grid ( igrid ), big_grid ( igrid ), &
pnode_grp, big_grid ( iref ) )
ENDIF
CALL dg_rho0_setup ( dg_coef ( igt ) % dg_rho0, small_grid ( igrid ) )
IF ( igrid /= igrid_old ) THEN
IF ( igrid /= iref ) THEN
IF ( pnode_grp % ionode ) THEN
WRITE ( iw, '( A,T63,3I6 )' ) &
' POL_GRIDS_INIT| grid number ', igrid
WRITE ( iw, '( A,T63,3I6 )' ) &
' POL_GRIDS_INIT| Small box grid ', small_grid ( igrid ) % npts
WRITE ( iw, '( A,T63,3I6 )' ) &
' POL_GRIDS_INIT| Full box grid ', big_grid ( igrid ) % npts
END IF
ELSE
IF ( pnode_grp % ionode ) THEN
WRITE ( iw, '( A,T63,3I6 )' ) &
' POL_GRIDS_INIT| Ref. grid number ', iref
WRITE ( iw, '( A,T63,3I6 )' ) &
' POL_GRIDS_INIT| Ref. small box grid ', small_grid ( iref ) % npts
WRITE ( iw, '( A,T63,3I6 )' ) &
' POL_GRIDS_INIT| Ref. full box grid ', big_grid ( iref ) % npts
END IF
ENDIF
ENDIF
igrid_old = igrid
END DO
! atom Loop
DO igt = 1, SIZE ( dg_atom )
ikind = dg_atom ( igt ) % dg_rho0 % kind
ki => rho0_basis_info ( ikind )
npts_s ( : ) = ki % number_of_grid_points
igrid = dg_atom ( igt ) % dg_rho0 % grid
cutoff_radius = dg_atom ( igt ) % dg_rho0 % cutoff_radius
CALL particle_list_grid ( part, dg_atom ( igt ) % plist, ki % atom_list, &
dg_atom ( igt ) % nparts, igrid )
IF ( igrid /= igrid_old ) THEN
IF ( igrid /= iref ) &
CALL dg_pol_grid_setup ( box, npts_s, cutoff_radius , &
small_grid ( igrid ), big_grid ( igrid ), &
pnode_grp, big_grid ( iref ) )
ENDIF
CALL dg_rho0_setup ( dg_atom ( igt ) % dg_rho0, small_grid ( igrid ) )
IF ( igrid /= igrid_old ) THEN
IF ( igrid /= iref ) THEN
IF ( pnode_grp % ionode ) THEN
WRITE ( iw, '( A,T63,3I6 )' ) &
' POL_GRIDS_INIT| grid number ', igrid
WRITE ( iw, '( A,T63,3I6 )' ) &
' POL_GRIDS_INIT| Small box grid ', small_grid ( igrid ) % npts
WRITE ( iw, '( A,T63,3I6 )' ) &
' POL_GRIDS_INIT| Full box grid ', big_grid ( igrid ) % npts
END IF
ELSE
IF ( pnode_grp % ionode ) THEN
WRITE ( iw, '( A,T63,3I6 )' ) &
' POL_GRIDS_INIT| Ref. grid number ', iref
WRITE ( iw, '( A,T63,3I6 )' ) &
' POL_GRIDS_INIT| Ref. small box grid ', small_grid ( iref ) % npts
WRITE ( iw, '( A,T63,3I6 )' ) &
' POL_GRIDS_INIT| Ref. full box grid ', big_grid ( iref ) % npts
END IF
ENDIF
ENDIF
igrid_old = igrid
END DO
END SUBROUTINE pol_grids_initialize
!******************************************************************************
SUBROUTINE assign_coefs_to_grids(ao,part,ki,ngt)
IMPLICIT NONE
TYPE (ao_type), INTENT (INOUT) :: ao
TYPE (particle_type), INTENT(IN), DIMENSION (:) :: part
TYPE (kind_info_type), DIMENSION (:), INTENT(IN) :: ki
INTEGER, INTENT(IN) :: ngt
! locals
INTEGER :: ngrid_types, grid_type_identifier, ikind, ii, ipart
INTEGER :: ngrid_kind, ncgf, icoef, icgf
ngrid_types = 0
ngrid_kind = 0
DO ikind = 1 , size ( ki )
ngrid_types = ngrid_types + ngrid_kind
DO ii = 1, ki ( ikind )%natom
ipart = ki ( ikind )% atom_list ( ii )
ncgf = ki ( ikind ) % orb_basis_set % ncgf
DO icgf = 1, ncgf
grid_type_identifier = ngrid_types + icgf
icoef = part ( ipart ) % coef_list( icgf )
ao % grid_id ( icoef ) = grid_type_identifier
ENDDO
END DO
ngrid_kind = ncgf
END DO
IF ( grid_type_identifier /= ngt .AND. ngt /= 0 ) THEN
CALL stop_program ( "assign_coef_to_grid", "inconsistent grid number" )
END IF
END SUBROUTINE assign_coefs_to_grids
!!*****
!******************************************************************************
END MODULE pol_grids
!******************************************************************************

349
src/pol_setup.F Normal file
View file

@ -0,0 +1,349 @@
!-----------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!***** cp2k/pol_setup [1.0] *
!!
!! NAME
!! pol_setup
!!
!! FUNCTION
!!
!! AUTHOR
!! gloria 080501
!!
!! MODIFICATION HISTORY
!!
!! NOTES
!!
!!
!! SOURCE
!******************************************************************************
MODULE pol_setup
USE atomic_kinds, ONLY : kind_info_type
USE basis_set_types, ONLY : read_basis_set, allocate_basis_set, &
deallocate_basis_set
USE force_fields, ONLY : ATOMNAMESLENGTH
USE global_types, ONLY: global_environment_type
USE kinds, ONLY : dbl
USE interactions, ONLY : exp_radius
USE molecule_types , ONLY : molecule_type, read_basis_type
USE particle_types, ONLY : particle_type
USE string_utilities, ONLY : str_search
USE termination, ONLY : stop_memory, stop_program
USE timings, ONLY : timeset, timestop
IMPLICIT NONE
PUBLIC :: read_pol_basis_sets, get_rcutsq_cgf
PRIVATE
CONTAINS
!******************************************************************************
SUBROUTINE read_pol_basis_sets(rho0_basis_info, drho_basis_info, &
mol_setup,part,atom_names,ncoef,globenv)
IMPLICIT NONE
TYPE ( kind_info_type ), POINTER, DIMENSION(:) :: rho0_basis_info
TYPE ( kind_info_type ), POINTER, DIMENSION(:) :: drho_basis_info
TYPE ( molecule_type ), DIMENSION ( : ), POINTER :: mol_setup
TYPE ( particle_type ), POINTER, DIMENSION(:) :: part
CHARACTER ( LEN = ATOMNAMESLENGTH ), DIMENSION ( : ), POINTER :: atom_names
integer, intent (out) :: ncoef
TYPE(global_environment_type), INTENT(IN) :: globenv
! TYPE (ao_type) :: ao
! locals
INTEGER :: i,nrho0,ndrho,ikind,ncgf,isos,imol, icount
TYPE ( read_basis_type ), POINTER, DIMENSION (:) :: pol_basis
TYPE ( read_basis_type ), POINTER, DIMENSION (:) :: rho0_basis
LOGICAL :: flag
! count the total number of basis functions for
! rho0 and polarization density
nrho0 = 0
ndrho = 0
icount = 0
DO imol=1, SIZE(mol_setup)
nrho0 = nrho0 + mol_setup(imol)%molpar%nbasis_rho0
ndrho = ndrho + mol_setup(imol)%molpar%nbasis_pol
END DO
ALLOCATE (pol_basis(ndrho), STAT = isos)
DO imol=1, SIZE(mol_setup)
DO i = 1, mol_setup(imol)%molpar%nbasis_pol
icount = icount + 1
pol_basis(icount)%aname = mol_setup(imol)%molpar%pol_basis(i)%aname
pol_basis(icount)%libname = mol_setup(imol)%molpar%pol_basis(i)%libname
pol_basis(icount)%eps = mol_setup(imol)%molpar%pol_basis(i)%eps
pol_basis(icount)%nmax = mol_setup(imol)%molpar%pol_basis(i)%nmax
END DO
END DO
ALLOCATE (rho0_basis(nrho0), STAT = isos)
icount = 0
DO imol=1, SIZE(mol_setup)
DO i = 1, mol_setup(imol)%molpar%nbasis_rho0
icount = icount + 1
rho0_basis(icount)%aname = mol_setup(imol)%molpar%rho0_basis(i)%aname
rho0_basis(icount)%libname = mol_setup(imol)%molpar%rho0_basis(i)%libname
rho0_basis(icount)%eps = mol_setup(imol)%molpar%rho0_basis(i)%eps
rho0_basis(icount)%nmax = mol_setup(imol)%molpar%rho0_basis(i)%nmax
END DO
END DO
! allocate the pointer drho_basis_info
NULLIFY (drho_basis_info)
CALL allocate_pol_basis_info(drho_basis_info,ndrho)
! allocate the pointer rho0_basis_info
NULLIFY (rho0_basis_info)
CALL allocate_pol_basis_info(rho0_basis_info,nrho0)
! fill the pointer drho_basis_info
CALL initialize_pol_basis_info(drho_basis_info,ndrho,pol_basis,part,atom_names, &
globenv)
! fill the pointer rho0_basis_info
CALL initialize_pol_basis_info(rho0_basis_info,nrho0,rho0_basis,part,atom_names, &
globenv)
! count the total number of basis functions (ncoef) needed to allocate
! the coefficents vectors
CALL get_number_of_coefs(drho_basis_info, part, ncoef)
END SUBROUTINE read_pol_basis_sets
! ***************************************************************************
SUBROUTINE allocate_pol_basis_info(ki,number_of_kinds)
INTEGER, INTENT(IN) :: number_of_kinds
TYPE (kind_info_type), POINTER, DIMENSION(:) :: ki
! *** Local variables ***
INTEGER :: ikind,istat,nkind
! ---------------------------------------------------------------------------
! *** Deallocate the old data structure for the atomic kGind information ***
IF (ASSOCIATED(ki)) THEN
nkind = size(ki)
DO ikind=1,nkind
IF (ASSOCIATED(ki(ikind)%aux_basis_set)) THEN
CALL deallocate_basis_set(ki(ikind)%aux_basis_set)
END IF
IF (ASSOCIATED(ki(ikind)%orb_basis_set)) THEN
CALL deallocate_basis_set(ki(ikind)%orb_basis_set)
END IF
DEALLOCATE (ki(ikind)%atom_list,STAT=istat)
IF (istat /= 0) CALL stop_memory("allocate_pol_basis_info", &
"ki(ikind)%atom_list")
DEALLOCATE (ki(ikind)%elec_conf,STAT=istat)
IF (istat /= 0) CALL stop_memory("allocate_pol_basis_info","ki(ikind)%elec_conf")
END DO
DEALLOCATE (ki,STAT=istat)
IF (istat /= 0) CALL stop_memory("allocate_pol_basis_info","ki")
END IF
nkind = number_of_kinds
IF (nkind >= 0) THEN
ALLOCATE (ki ( nkind ),STAT=istat)
IF (istat /= 0) CALL stop_memory("allocate_pol_basis_info","ki(nkind)",0)
DO ikind=1,nkind
NULLIFY ( ki ( ikind ) % aux_basis_set )
NULLIFY ( ki ( ikind ) % orb_basis_set )
NULLIFY ( ki ( ikind ) % atom_list )
NULLIFY ( ki ( ikind ) % elec_conf )
END DO
ELSE
CALL stop_program("allocate_pol_basis_info",&
"invalid number of atomic kinds")
END IF
END SUBROUTINE allocate_pol_basis_info
! *****************************************************************************
SUBROUTINE initialize_pol_basis_info(ki,basis_size,molinfo,part,atom_names,globenv)
IMPLICIT NONE
TYPE(kind_info_type), POINTER, DIMENSION(:) :: ki
INTEGER :: basis_size
TYPE(particle_type), POINTER, DIMENSION(:) :: part
TYPE(read_basis_type), DIMENSION(:), POINTER :: molinfo
CHARACTER ( LEN = ATOMNAMESLENGTH ), DIMENSION ( : ), POINTER :: atom_names
TYPE(global_environment_type), INTENT(IN) :: globenv
! locals
CHARACTER(len=20) :: string
INTEGER :: icount, itype, nt, ikind
nt = size (atom_names)
DO ikind= 1, size(ki)
ki(ikind)%element_symbol = molinfo(ikind)%aname
ki(ikind)%orb_basis_set_name = molinfo(ikind)%libname
ki(ikind)%number_of_grid_points = molinfo(ikind)%nmax
string=molinfo(ikind)%aname
itype = str_search(atom_names,nt,string)
CALL get_atom_list(ki(ikind),part,itype)
CALL allocate_basis_set(ki(ikind)%orb_basis_set)
CALL read_basis_set(ki(ikind)%element_symbol,ki(ikind)%orb_basis_set_name, &
ki(ikind)%orb_basis_set, globenv)
CALL get_radii(ki(ikind), moliNFo(ikind)%eps)
ENDDO
END SUBROUTINE initialize_pol_basis_info
!------------------------------------------------------------------------------!
SUBROUTINE get_atom_list(ki,part,itype)
IMPLICIT NONE
TYPE(kind_info_type), INTENT(inout) :: ki
TYPE(particle_type), DIMENSION(:), INTENT(inout) :: part
INTEGER, INTENT(in) :: itype
! locals
INTEGER :: isos,i,icount
icount=0
DO i=1, SIZE(part)
IF (part(i)%prop%ptype==itype) icount=icount+1
ENDDO
ki%natom = icount
ALLOCATE(ki%atom_list(icount), STAT=isos)
IF (isos/=0) CALL stop_memory('get_atom_list', 'atom_list')
icount=0
DO i=1, SIZE(part)
IF (part(i)%prop%ptype==itype) THEN
icount=icount+1
ki%atom_list(icount)=i
ENDIF
ENDDO
END SUBROUTINE get_atom_list
!------------------------------------------------------------------------------!
SUBROUTINE get_radii(ki,eps)
IMPLICIT NONE
TYPE (kind_info_type), INTENT(inout) :: ki
REAL (dbl), INTENT(inout) :: eps
!locals
REAL(dbl) :: gcca,kind_radius,pgf_radius,set_radius,zeta
INTEGER :: ipgf,ikind,iset,ishell,l
!---------------------------------------------------------------------------
kind_radius = 0.0_dbl
DO iset=1,ki%orb_basis_set%nset
set_radius = 0.0_dbl
DO ipgf=1,ki%orb_basis_set%npgf(iset)
pgf_radius = 0.0_dbl
DO ishell=1,ki%orb_basis_set%nshell(iset)
l = ki%orb_basis_set%l(ishell,iset)
gcca = ki%orb_basis_set%gcc(ipgf,ishell,iset)
zeta = ki%orb_basis_set%zet(ipgf,iset)
pgf_radius = MAX(pgf_radius,exp_radius(l,zeta,eps,gcca))
END DO
ki%orb_basis_set%pgf_radius(ipgf,iset) = pgf_radius
set_radius = MAX(set_radius,pgf_radius)
END DO
ki%orb_basis_set%set_radius(iset) = set_radius
kind_radius = MAX(kind_radius,set_radius)
END DO
ki%orb_basis_set%kind_radius = kind_radius
END SUBROUTINE get_radii
!------------------------------------------------------------------------------
SUBROUTINE get_number_of_coefs(ki, part, ncoefs)
IMPLICIT NONE
TYPE (kind_info_type), INTENT(in), DIMENSION (:) :: ki
TYPE (particle_type), INTENT(inout), DIMENSION (:) :: part
INTEGER, INTENT (OUT) :: ncoefs
! locals
INTEGER :: icount, ii, ikind, i, icgf, ncgf, isos
ncoefs = 0
icount = 0
DO ikind = 1, size(ki)
DO ii = 1, ki(ikind)%natom
i = ki(ikind)%atom_list(ii)
NULLIFY (part(i) % coef_list)
ncgf = ki(ikind)% orb_basis_set % ncgf
ncoefs = ncoefs + ki(ikind)% orb_basis_set % ncgf
ALLOCATE(part (i) %coef_list(ncgf), STAT=isos)
IF (isos/=0) CALL stop_memory('get_number_of_coef', 'coeflist')
DO icgf = 1, ncgf
icount = icount + 1
part (i) % coef_list(icgf)=icount
ENDDO
END DO
END DO
END SUBROUTINE get_number_of_coefs
!------------------------------------------------------------------------------
SUBROUTINE get_rcutsq_cgf(ki,rsq)
IMPLICIT NONE
TYPE (kind_info_type), DIMENSION (:), INTENT(in) :: kI
REAL (dbl), INTENT(out), DIMENSION (:,:) :: rsq
!locals
INTEGER :: i, j
REAL (dbl) :: rcut
! calculate the cutoff radius for the Gaussians needed for the overlap list
DO i = 1, size(ki)
DO j = 1, size (ki)
rcut = ki(i)%orb_basis_set%kind_radius + &
ki(j)%orb_basis_set%kind_radius
rsq (i,j) = rcut * rcut
END DO
END DO
END SUBROUTINE get_rcutsq_cgf
!------------------------------------------------------------------------------!
END MODULE pol_setup
!------------------------------------------------------------------------------!

2
src/pw_grids.F
View file

@ -405,7 +405,7 @@ SUBROUTINE pw_grid_setup ( cell, pw_grid, cutoff, pe_group, info, fft_usage, &
IF ( pw_grid % grid_span /= ref_grid % grid_span ) THEN
CALL stop_program ( "pw_grid_setup", "Incompatible grid types" )
END IF
IF ( pw_grid % spherical .EQV. ref_grid % spherical ) THEN
IF ( pw_grid % spherical .NEQV. ref_grid % spherical ) THEN
CALL stop_program ( "pw_grid_setup", "Incompatible cutoff schemes" )
END IF
END IF

71
src/pw_types.F
View file

@ -2,7 +2,7 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****s* cp2k/pw_types [1.0] *
!!***** cp2k/pw_types [1.0] *
!!
!! NAME
!! pw_types
@ -593,7 +593,8 @@ END SUBROUTINE pw_transfer
!! apsi
!!
!! MODIFICATION HISTORY
!! none
!! CJM (23-May-01) : check for pw_grid % identifier, allow copy if
!! in_use == COMPLEXDATA1D and in_space == RECIPROCALSPACE
!!
!! NOTES
!! Currently only adding of respective types allowed,
@ -611,33 +612,59 @@ SUBROUTINE pw_add ( pw1, pw2, pw3 )
TYPE ( pw_type ), INTENT ( INOUT ) :: pw3
! Local
INTEGER :: handle
INTEGER :: handle, ng3, ng2, ng1, ng
REAL ( dbl ) :: flop
!------------------------------------------------------------------------------
CALL timeset ( 'PW_ADD', 'I', 'Mflops', handle )
IF ( pw3 % pw_grid % identifier /= pw2 % pw_grid % identifier .OR. &
pw3 % pw_grid % identifier /= pw1 % pw_grid % identifier ) THEN
IF ( pw1 % in_use == REALDATA1D .AND. pw2 % in_use == REALDATA1D .AND. &
pw3 % in_use == REALDATA1D ) THEN
pw3 % cr = pw1 % cr + pw2 % cr
flop = REAL ( SIZE ( pw2 % cr ), dbl )
ELSE IF ( pw1 % in_use == COMPLEXDATA1D .AND. &
pw2 % in_use == COMPLEXDATA1D .AND. pw3 % in_use == COMPLEXDATA1D ) THEN
pw3 % cc = pw1 % cc + pw2 % cc
flop = REAL ( 2 * SIZE ( pw2 % cc ), dbl )
ELSE IF ( pw1 % in_use == REALDATA3D .AND. pw2 % in_use == REALDATA3D .AND. &
pw3 % in_use == REALDATA3D ) THEN
pw3 % cr3d = pw1 % cr3d + pw2 % cr3d
flop = REAL ( SIZE ( pw2 % cr3d ), dbl )
ELSE IF ( pw1 % in_use == COMPLEXDATA3D .AND. &
pw2 % in_use == COMPLEXDATA3D .AND. pw3 % in_use == COMPLEXDATA3D ) THEN
pw3 % cc3d = pw1 % cc3d + pw2 % cc3d
flop = REAL ( 2 * SIZE ( pw2 % cc3d ), dbl )
ELSE
CALL stop_program ( "pw_add", "no suitable data field" )
IF ( pw1 % in_use == COMPLEXDATA1D .AND. &
pw2 % in_use == COMPLEXDATA1D .AND. &
pw3 % in_use == COMPLEXDATA1D .AND. &
pw1 % in_space == RECIPROCALSPACE .AND. &
pw2 % in_space == RECIPROCALSPACE ) THEN
ng1 = SIZE ( pw1 % cc )
ng2 = SIZE ( pw2 % cc )
ng3 = SIZE ( pw3 % cc )
IF ( ( ng3 < ng2 ) .OR. ( ng3 < ng1 ) ) &
CALL stop_program ( "pw_add", &
"final data field smaller than component data fields" )
ng = MIN ( ng1, ng2 )
pw3 % cc ( 1:ng ) = pw2 % cc ( 1:ng ) + pw1 % cc ( 1:ng )
IF ( ng1 > ng ) pw3 % cc ( ng+1:ng1 ) = pw1 % cc ( ng+1:ng1 )
IF ( ng2 > ng ) pw3 % cc ( ng+1:ng2 ) = pw2 % cc ( ng+1:ng2 )
IF ( ( ng3 > ng1 ) .AND. ( ng3 > ng2 ) ) THEN
ng = MAX ( ng1, ng2 )
pw3 % cc ( ng+1:ng3 ) = CMPLX ( 0._dbl, 0._dbl, dbl )
ENDIF
flop = ng3
ELSE
CALL stop_program ( "pw_add", "no suitable data field" )
END IF
ELSE
IF ( pw1 % in_use == REALDATA1D .AND. pw2 % in_use == REALDATA1D .AND. &
pw3 % in_use == REALDATA1D ) THEN
pw3 % cr = pw1 % cr + pw2 % cr
flop = REAL ( SIZE ( pw2 % cr ), dbl )
ELSE IF ( pw1 % in_use == COMPLEXDATA1D .AND. &
pw2 % in_use == COMPLEXDATA1D .AND. pw3 % in_use == COMPLEXDATA1D ) THEN
pw3 % cc = pw1 % cc + pw2 % cc
flop = REAL ( 2 * SIZE ( pw2 % cc ), dbl )
ELSE IF ( pw1 % in_use == REALDATA3D .AND. pw2 % in_use == REALDATA3D .AND. &
pw3 % in_use == REALDATA3D ) THEN
pw3 % cr3d = pw1 % cr3d + pw2 % cr3d
flop = REAL ( SIZE ( pw2 % cr3d ), dbl )
ELSE IF ( pw1 % in_use == COMPLEXDATA3D .AND. &
pw2 % in_use == COMPLEXDATA3D .AND. pw3 % in_use == COMPLEXDATA3D ) THEN
pw3 % cc3d = pw1 % cc3d + pw2 % cc3d
flop = REAL ( 2 * SIZE ( pw2 % cc3d ), dbl )
ELSE
CALL stop_program ( "pw_add", "no suitable data field" )
END IF
END IF
flop = flop * 1.e-6_dbl
CALL timestop ( flop, handle )

4
src/pw_xc.F
View file

@ -32,9 +32,9 @@ MODULE pw_xc
NOSPACE, REALSPACE, RECIPROCALSPACE
USE functionals, ONLY : xc_type, becke88, lyp, pade, pbe_c, pbe_x, &
perdew86_x, perdew86_c, slater, vwn_c, vwn_x
IMPLICIT NONE
INTERFACE pw_calculate_xc
MODULE PROCEDURE pw_calculate_xcns, pw_calculate_xcspin
END INTERFACE

48
src/pws.F
View file

@ -110,9 +110,9 @@ SUBROUTINE pw_poisson_solver ( density, ehartree, vhartree, dvhartree, h_stress
IMPLICIT NONE
! Arguments
TYPE ( coeff_type ), INTENT ( INOUT ), TARGET :: density
REAL ( dbl ) , INTENT ( OUT ) :: ehartree
TYPE ( coeff_type ), INTENT ( OUT ) :: vhartree
TYPE ( coeff_type ), INTENT ( IN ) :: density
REAL ( dbl ) , INTENT ( OUT ), OPTIONAL :: ehartree
TYPE ( coeff_type ), INTENT ( OUT ), TARGET, OPTIONAL :: vhartree
TYPE ( coeff_type ), DIMENSION ( 3 ), INTENT ( OUT ), OPTIONAL :: dvhartree
REAL ( dbl ), DIMENSION ( 3, 3 ), INTENT ( OUT ), OPTIONAL :: h_stress
@ -121,35 +121,40 @@ SUBROUTINE pw_poisson_solver ( density, ehartree, vhartree, dvhartree, h_stress
TYPE ( pw_type ) :: rhog, dvg ( 3 )
REAL ( dbl ) :: ffa
INTEGER, DIMENSION ( 3 ) :: n
INTEGER :: i, gpt, alpha, beta
INTEGER :: i, gpt, alpha, beta, ng
!------------------------------------------------------------------------------
pw_grid => density % pw % pw_grid
IF ( .NOT. ASSOCIATED ( gf ) ) THEN
CALL stop_program ( "pw_poisson_solver", "init_pw_poisson_solver "// &
"has to be used before pw_poisson_solver" )
END IF
! density in G space
CALL pw_allocate ( rhog, pw_grid, COMPLEXDATA1D )
rhog % in_space = RECIPROCALSPACE
! update the greens function
CALL pw_green_fn ( gf )
IF ( PRESENT ( vhartree ) ) THEN
pw_grid => vhartree % pw % pw_grid
! density in G space
CALL pw_allocate ( rhog, pw_grid, COMPLEXDATA1D )
rhog % in_space = RECIPROCALSPACE
SELECT CASE ( gf % method )
CASE ( PERIODIC3D, ANALYTIC2D, ANALYTIC1D, ANALYTIC0D, MT2D, MT1D, MT0D )
CALL pw_transfer ( density % pw, rhog )
rhog % cc = rhog % cc * gf % influence_function % cr
CALL pw_transfer ( rhog, vhartree % pw )
CASE ( HOCKNEY2D, HOCKNEY1D, HOCKNEY0D )
CALL stop_program ( "pw_poisson_solver", " Method not implemented ")
CASE DEFAULT
CALL stop_program ( "pw_poisson_solver", " Method not implemented ")
END SELECT
! update the greens function
CALL pw_green_fn ( gf )
ng = SIZE ( pw_grid % gsq )
SELECT CASE ( gf % method )
CASE ( PERIODIC3D, ANALYTIC2D, ANALYTIC1D, ANALYTIC0D, MT2D, MT1D, MT0D )
CALL pw_transfer ( density % pw, rhog )
rhog % cc ( 1 : ng ) = rhog % cc ( 1 : ng ) * &
gf % influence_function % cr ( 1 : ng )
CALL pw_transfer ( rhog, vhartree % pw )
CASE ( HOCKNEY2D, HOCKNEY1D, HOCKNEY0D )
CALL stop_program ( "pw_poisson_solver", " Method not implemented ")
CASE DEFAULT
CALL stop_program ( "pw_poisson_solver", " Method not implemented ")
END SELECT
ENDIF
! do we need to calculate the energy
IF ( PRESENT ( ehartree ) ) &
ehartree = 0.5_dbl * pw_integral_ab ( density % pw, vhartree % pw )
! do we need to calculate the derivative of the potential?
@ -193,6 +198,7 @@ SUBROUTINE pw_poisson_solver ( density, ehartree, vhartree, dvhartree, h_stress
END DO
END IF
IF ( PRESENT ( vhartree ) ) &
CALL pw_deallocate ( rhog )
END SUBROUTINE pw_poisson_solver

2
src/radial_solver.F
View file

@ -96,7 +96,7 @@ CONTAINS
ab ( 4, i ) = ab ( 4, i ) * ww
END DO
CALL dgbsv ( n, 1, 1, 1, ab, 4, ipiv, s, n, info )
CALL DGBSV ( n, 1, 1, 1, ab, 4, ipiv, s, n, info )
IF ( info /= 0 ) call stop_program ( "numerov", "DGBSV: info /= 0" )
g ( 1:n ) = s ( 1:n )

10
src/spme.F
View file

@ -26,7 +26,7 @@ MODULE spme
PW_REALDATA3D, PW_COMPLEXDATA1D, &
PW_REALSPACE, PW_RECIPROCALSPACE
USE dg_types, ONLY : dg_type
USE dgs, ONLY : dg_sum_patch, dg_sum_patch_force
USE dgs, ONLY : dg_sum_patch, dg_sum_patch_force_3d
USE ewald_parameters_types, ONLY : ewald_parameters_type
USE greens_fn, ONLY : greens_function_type
USE hartree, ONLY : calculate_hartree
@ -215,7 +215,7 @@ SUBROUTINE spme_evaluate ( dg, part, fg_coulomb, vg_coulomb, pv_g, box, &
CALL get_patch ( part, box, grid_spme % npts, p1, rhos )
! add boxes to real space grid (big box)
CALL dg_sum_patch_force ( drpot, rhos, centre(:,p1), fat )
CALL dg_sum_patch_force_3d ( drpot, rhos, centre(:,p1), fat )
fg_coulomb ( 1, p1 ) = fg_coulomb ( 1, p1 ) - fat ( 1 ) * dvols
fg_coulomb ( 2, p1 ) = fg_coulomb ( 2, p1 ) - fat ( 2 ) * dvols
fg_coulomb ( 3, p1 ) = fg_coulomb ( 3, p1 ) - fat ( 3 ) * dvols
@ -609,7 +609,7 @@ SUBROUTINE spme_coeff_calculate ( n, coeff )
! Locals
INTEGER :: i, j, l, m
REAL ( dbl ), DIMENSION ( n, -(n-1):n-1, 0:n-1 ) :: a
REAL ( dbl ), DIMENSION ( n, -n:n, 0:n-1 ) :: a
REAL ( dbl ) :: b
!------------------------------------------------------------------------------
@ -639,10 +639,6 @@ SUBROUTINE spme_coeff_calculate ( n, coeff )
coeff ( j, i ) = a ( n, j, i )
END DO
END DO
!!!!!!!!!
!This seems to cause the SGI compiler to work correctly!!!!!
write(*,*)
!!!!!!!!!
END SUBROUTINE spme_coeff_calculate

12
src/structure_types.F
View file

@ -14,6 +14,7 @@
!!
!! MODIFICATION HISTORY
!! Harald Forbert (Dec-2000): Add ll-data to structure_type
!! CJM 14-April-2001 : added coeff_type to structure_type
!!
!! SOURCE
!******************************************************************************
@ -21,6 +22,8 @@
MODULE structure_types
USE kinds, ONLY : dbl
USE atomic_kinds, ONLY : kind_info_type
USE coefficient_types, ONLY : coeff_type
USE molecule_types, ONLY : molecule_structure_type, particle_node_type
USE particle_types, ONLY : particle_type
USE simulation_cell, ONLY : cell_type
@ -46,7 +49,7 @@ MODULE structure_types
REAL ( dbl ) :: ismeff
REAL ( dbl ), DIMENSION ( 2 , 2 ) :: harmonic_prop
END TYPE pimd_data_type
TYPE structure_type
CHARACTER ( LEN = 30 ) :: name
TYPE ( cell_type ) :: box
@ -54,8 +57,13 @@ MODULE structure_types
TYPE ( particle_type ), POINTER :: part ( : )
TYPE ( particle_node_type ), POINTER :: pnode ( : )
TYPE ( molecule_structure_type ), POINTER :: molecule ( : )
TYPE ( linklist_internal_data_type) :: ll_data
TYPE ( linklist_internal_data_type), POINTER :: ll_data ( : )
TYPE ( pimd_data_type ) :: pimd_data
TYPE ( kind_info_type ), POINTER :: rho0_basis_info ( : )
TYPE ( kind_info_type ), POINTER :: drho_basis_info ( : )
TYPE ( coeff_type ), POINTER, DIMENSION ( : ) :: coef_pos
TYPE ( coeff_type ), POINTER, DIMENSION ( : ) :: coef_vel
TYPE ( coeff_type ), POINTER, DIMENSION ( : ) :: coef_force
END TYPE structure_type
TYPE interaction_type

123
src/timings.F
View file

@ -28,20 +28,19 @@ MODULE timings
USE termination, ONLY : stop_program
USE timesl, ONLY : cputime, walltime
USE util, ONLY : sort
USE message_passing, ONLY : mp_world_setup
PRIVATE
PUBLIC :: timeset, timestop, timeprint, trace_debug, print_stack
! subroutine calling stack
!!subroutine calling stack
INTEGER, PARAMETER :: max_stack = 100
INTEGER :: stack_size
CHARACTER ( LEN = 31 ) :: routine_stack ( max_stack ), routine_name
CHARACTER ( LEN = 30 ) :: routine_stack ( max_stack ), routine_name
! subroutine timing and performance information
!!subroutine timing and performance information
INTEGER, PARAMETER :: max_timer = 200
CHARACTER ( LEN = 31 ) :: clock_name ( max_timer )
CHARACTER ( LEN = 30 ) :: clock_name ( max_timer )
INTEGER :: routine_calls ( max_timer )
INTEGER :: clock_status ( max_timer )
REAL ( dbl ) :: clock_start ( max_timer )
@ -55,7 +54,6 @@ MODULE timings
INTEGER :: iod = 0
LOGICAL :: trace = .false.
CHARACTER ( LEN = 15 ) :: astring, bstring
!!*****
!******************************************************************************
@ -101,7 +99,7 @@ SUBROUTINE timeset ( name, ct, pu, handle )
!..if debug mode echo the subroutine name
IF ( trace ) write ( iod, '(A,A,A)') &
astring,TRIM(name)," started ..."
" >>>>> ",TRIM(name)," started ..."
!..is this a new timer or has this routine been called before
handle = 0
@ -188,7 +186,7 @@ SUBROUTINE timestop ( perf, handle )
!..if debug mode echo the subroutine name
IF ( trace ) write ( iod, '(A,A,A)') &
bstring,TRIM(routine_stack(stack_size))," ... ended"
" <<<<< ",TRIM(routine_stack(stack_size))," ... ended"
!..remove name from subroutine stack
stack_size = stack_size - 1
@ -215,74 +213,84 @@ END SUBROUTINE timestop
!! SOURCE
!******************************************************************************
SUBROUTINE timeprint(iw,global_print_level)
SUBROUTINE timeprint ( iw, global_print_level )
IMPLICIT NONE
! Arguments
INTEGER, INTENT(IN) :: global_print_level,iw
INTEGER :: iw, global_print_level
! Locals
REAL(dbl) :: maxtime,mintime,perf
INTEGER :: i,j
CHARACTER(LEN=1), DIMENSION(0:1) :: ct = (/"I","E"/)
INTEGER, DIMENSION(max_timer) :: index
INTEGER :: index ( max_timer ), i, j, ia, ie
REAL ( dbl ) :: maxtime, mintime, perf
CHARACTER ( LEN = 17 ) :: sname
CHARACTER ( LEN = 6 ) :: punit
CHARACTER ( LEN = 1 ) :: blank = ' ', ct ( 0 : 1 ) = (/ 'I', 'E' /)
!------------------------------------------------------------------------------
IF ((num_timer > 0).AND.(global_print_level >= 0)) THEN
IF ( num_timer > 0 .AND. global_print_level >= 0 ) THEN
CALL sort(clock_accu,num_timer,index)
CALL sort ( clock_accu, num_timer, index )
maxtime = clock_accu(num_timer)
maxtime = clock_accu ( index ( num_timer ) )
IF (global_print_level > 10) THEN
IF ( global_print_level > 10 ) THEN
mintime = zero
ELSE IF (global_print_level > 5) THEN
mintime = maxtime*0.002_dbl
ELSE IF ( global_print_level > 5 ) THEN
mintime = maxtime * 0.002_dbl
ELSE
mintime = maxtime*0.02_dbl
mintime = maxtime * 0.02_dbl
END IF
WRITE (UNIT=iw,FMT="(/,T2,A)") REPEAT("-",79)
WRITE (UNIT=iw,FMT="(T2,A,T80,A)") "-","-"
WRITE (UNIT=iw,FMT="(T2,A,T35,A,T80,A)") "-","T I M I N G","-"
WRITE (UNIT=iw,FMT="(T2,A,T80,A)") "-","-"
WRITE (UNIT=iw,FMT="(T2,A)") REPEAT("-",79)
WRITE (UNIT=iw,FMT="(T2,A,T36,A,T43,A,T55,A,T69,A)")&
"SUBROUTINE","CALLS","CPU TIME T","PERFORMANCE","ELAPSED TIME"
WRITE ( iw, '( /, 1X, 79("-") )' )
WRITE ( iw, '( " -", 77X, "-" )' )
WRITE ( iw, '( " -", 32X, A, 32X, "-" )' ) ' T I M I N G '
WRITE ( iw, '( " -", 77X, "-" )' )
WRITE ( iw, '( 1X, 79("-"), / )' )
WRITE ( iw, '( A, A, A )' ) ' SUBROUTINE', ' CALLS ', &
' CPU TIME T PERFORMANCE ELAPSED TIME'
DO i=num_timer,1,-1
DO i = num_timer, 1, -1
j = index(i)
j = index ( i )
IF ( clock_accu ( j ) > mintime ) THEN
IF (clock_accu(i) > mintime) THEN
IF ( clock_accu ( i ) > zero ) THEN
perf = perf_count ( j ) / clock_accu ( i )
ELSE
perf = zero
END IF
IF (clock_accu(i) > zero) THEN
perf = perf_count(j)/clock_accu(i)
ELSE
perf = zero
END IF
sname = blank
ia = 0
ie = 0
CALL xstring ( clock_name ( j ), ia, ie )
ie = min ( ia+16-1, ie )
sname ( 1 : ie-ia+1 ) = clock_name ( j ) ( ia : ie )
punit = blank
IF (LEN_TRIM(perf_unit(j)) == 0) THEN
WRITE (UNIT=iw,&
FMT="(T2,A,T34,I7,T42,F9.2,T52,A,T60,A,T67,F13.2)")&
TRIM(clock_name(j)),routine_calls(j),clock_accu(i),&
ct(clock_type(j)),"---",wallclock_accu(j)
ELSE
WRITE (UNIT=iw,&
FMT="(T2,A,T34,I7,T42,F9.2,T52,A,T54,F5.1,T60,A,T67,F13.2)")&
TRIM(clock_name(j)),routine_calls(j),clock_accu(i),&
ct(clock_type(j)),perf,TRIM(perf_unit(j)),wallclock_accu(j)
END IF
END IF
CALL xstring ( perf_unit ( j ), ia, ie )
IF ( ia > ie ) THEN
punit = '--- '
WRITE ( iw, &
'( 1X, A16, I7, F15.2, 1X, A1, 12X, 1X, A6, 5X, F15.2 )' ) &
sname, routine_calls ( j ), clock_accu ( i ), &
ct ( clock_type ( j ) ), punit, wallclock_accu ( j )
ELSE
ie = MIN ( ia+6-1, ie )
punit ( 6-ie+ia : 6 ) = perf_unit ( j ) ( ia : ie )
WRITE ( iw, &
'( 1X, A16, I7, F15.2, 1X, A1, F12.2, 1X, A6, 5X, F15.2 )' &
) &
sname, routine_calls ( j ), clock_accu ( i ), &
ct ( clock_type ( j ) ), perf, punit, wallclock_accu ( j )
END IF
END IF
END DO
WRITE (UNIT=iw,FMT="(T2,A,/)") REPEAT("-",79)
WRITE ( iw, '( 1X, 79("-"), / )' )
END IF
END SUBROUTINE timeprint
@ -316,7 +324,6 @@ SUBROUTINE trace_debug ( set, iw )
! Locals
INTEGER :: ia, ie
INTEGER :: numtask, taskid, groupid
CALL uppercase(set)
CALL xstring ( set, ia, ie )
@ -328,16 +335,6 @@ SUBROUTINE trace_debug ( set, iw )
call stop_program ("trace_debug","unknown action")
ENDIF
IF ( PRESENT ( iw ) ) iod = iw
#if defined (__parallel)
CALL mp_world_setup ( numtask, taskid, groupid )
astring = " >>>>>>>>>>>>> "
bstring = " >>>>>>>>>>>>> "
WRITE(astring(1:9),'(i4,a,i4)') taskid,":",numtask
WRITE(bstring(1:9),'(i4,a,i4)') taskid,":",numtask
#else
astring = " >>>>>>>>>>>>> "
bstring = " >>>>>>>>>>>>> "
#endif
END SUBROUTINE trace_debug

1
src/util.F
View file

@ -56,7 +56,6 @@ SUBROUTINE sort2 ( arr, n, index )
REAL ( dbl ) :: a, temp
!------------------------------------------------------------------------------
DO i = 1, n
index(i) = i
END DO