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

update to version 01/01/2001

Browse source 
svn-origin-rev: 11
This commit is contained in:
Jürg Hutter 2001-06-26 18:06:45 +00:00
parent e221094355
commit e7b68385e6
127 changed files with 5834 additions and 2178 deletions
Download patch file Download diff file Expand all files Collapse all files

2
src/MACHINEDEFS.DEC
View file

@ -8,7 +8,7 @@ FCfixed = f95 -fixed
LDR = f95
AR = ar -r
CDEFS = -D__DEC -D__FFTW
CDEFS = -D__DEC -D__FFTW -D__FFTSG
DEBUG = -g -C
OPT = -w0 -O5 -fast -ansi_alias -ansi_args -fp_reorder -tune host -std1

2
src/MACHINEDEFS.IBM
View file

@ -8,7 +8,7 @@ FCfixed = xlf90 -qfixed
LDR = xlf90
AR = ar -r
CDEFS = -D__AIX -D__FFTW
CDEFS = -D__AIX -D__FFTW -D__FFTSG -D__FFTESSL
DEBUG = -g -C -d -qextchk -qinitauto=AA -qlanglvl=90pure \
-qflttrap=overflow:zerodivide:invalid:enable \

2
src/MACHINEDEFS.PGI
View file

@ -8,7 +8,7 @@ FCfixed = pgf90 -Mfixed
LDR = pgf90
AR = ar -r
CDEFS = -D__PGI -D__FFTW
CDEFS = -D__PGI -D__FFTW -D__FFTSG
DEBUG = -g -C
OPT = -O2 -fast

2
src/MACHINEDEFS.T3E
View file

@ -8,7 +8,7 @@ FCfixed = f90 -f fixed
LDR = f90
AR = ar -r
CDEFS = -D__T3E -D__FFTW
CDEFS = -D__T3E -D__FFTW -D__FFTSG
DEBUG = -g -R abcps -ei -m 0
OPT = -O3 -eIn -m 1

2
src/Makefile
View file

@ -1,6 +1,6 @@
.SUFFIXES: .o .F .d
CP2KHOME= ${HOME}/CP2K/
CP2KHOME= ${HOME}/CP2K
FORPAR = $(CP2KHOME)/tools/forpar.x -chkint
SFMAKEDEPEND = $(CP2KHOME)/tools/sfmakedepend -m int -s -f

14
src/OBJECTDEFS
View file

@ -1,29 +1,31 @@
OBJECTS_GENERIC = \
atoms_input.o band.o brillouin.o cntl_input.o \
atom_module.o atoms_input.o band.o brillouin.o cntl_input.o \
coefficients.o coefficient_types.o \
constraint.o control_module.o convert_units.o \
dgs.o dg_rho0s.o dg_types.o \
dump.o cp2k.o cp2k_input.o environment.o \
eigenvalueproblems.o erf_fn.o global_types.o header.o \
ewalds.o ewald_parameters_types.o \
fermi.o fft_tools.o fftsg_lib.o fftw_lib.o \
fist_debug.o fist_global.o \
fermi.o fft_tools.o fftessl_lib.o fftsg_lib.o fftw_lib.o \
fist_debug.o \
fist_force.o fist_force_numer.o \
fist_intra_force.o fist_nonbond_force.o force_control.o \
force_fields.o integrator.o input_types.o \
initialize_particle_types.o initialize_molecule_types.o \
initialize_extended_types.o io_parameters.o lapack.o \
initialize_extended_types.o lapack.o \
k290.o kinds.o kpoint_initialization.o kpoints.o \
library_tests.o linklists.o linklist_control.o \
linklist_cell_types.o linklist_cell_list.o \
linklist_types.o \
linklist_utilities.o linklist_verlet_list.o machine.o \
mathconstants.o mol_force.o molecule_types.o \
molecule_input.o md.o message_passing.o nhc.o nose.o \
pair_potential.o particle_lists.o particle_types.o \
periodic_table.o parallel.o parallel_include.o parser.o \
physcon.o pme.o pws.o pw_types.o pw_grids.o \
pw_grid_types.o simulation_cell.o splines.o \
pw_grid_types.o radial_grids.o radial_poisson.o radial_solver.o \
radial_operators.o radial_util.o simulation_cell.o splines.o \
stop_program.o string_utilities.o structure_factors.o \
structure_factor_types.o structure_types.o timesl.o \
tbmd_debug.o tbmd_force.o tbmd_global.o tbmd_initialize.o \
tbmd_debug.o tbmd_force.o tbmd_initialize.o \
tbmd_input.o tbmd_types.o timings.o unit.o util.o

31
src/amoeba.F
View file

@ -2,6 +2,20 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 1999 MPI fuer Festkoerperforschung, Stuttgart !
!-----------------------------------------------------------------------------!
!!****** cp2k/amoeba [1.0] *
!!
!! NAME
!! amoeba, NUMERICAL RECIPIES SUBROUTINE AMOEBA ADAPTED FOR USE IN CP2K
!!
!! FUNCTION
!!
!! AUTHOR
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE amoeba
@ -13,8 +27,24 @@ MODULE amoeba
PRIVATE
PUBLIC :: amoeba_evaluate
!!*****
CONTAINS
!!****** amoeba/amoeba_evaluate [1.0] *
!!
!! NAME
!! amoeba_evaluate
!!
!! FUNCTION
!!
!! AUTHOR
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
!-----------------------------------------------------------------------------!
! NUMERICAL RECIPIES SUBROUTINE AMOEBA ADAPTED FOR USE IN CP2K !
!-----------------------------------------------------------------------------!
@ -115,6 +145,7 @@ CONTAINS
END SUBROUTINE amoeba_evaluate
!!*****
!******************************************************************************
END MODULE amoeba

156
src/atom_module.F Normal file
View file

@ -0,0 +1,156 @@
!-----------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/atom_module [1.0] *
!!
!! NAME
!! atom_module
!!
!! FUNCTION
!!
!! AUTHOR
!! JGH, 17-NOV-2000
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE atom_module
USE global_types, ONLY : global_environment_type
USE header, ONLY : atom_header
USE kinds, ONLY : dbl
USE mathconstants, ONLY : zero, pi
USE timings, ONLY : timeset, timestop
!deb
USE radial_grids, ONLY : radial_grid_type, init_radial_grid
USE radial_operators, ONLY : rad_kinetic
USE radial_poisson, ONLY : rg_poisson
USE radial_util, ONLY : radial_integrate, radial_diff
!deb
IMPLICIT NONE
PRIVATE
PUBLIC :: atomic_code
!!*****
!-----------------------------------------------------------------------------!
CONTAINS
!-----------------------------------------------------------------------------!
!!****** atom_module/atomic_code [1.0] *
!!
!! NAME
!! atomic_code
!!
!! FUNCTION
!!
!! AUTHOR
!! JGH, 17-NOV-2000
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
SUBROUTINE atomic_code ( globenv )
IMPLICIT NONE
!Arguments
TYPE ( global_environment_type ), INTENT ( INOUT ) :: globenv
!Locals
REAL ( dbl ) :: xx
INTEGER :: handle
!deb
REAL ( dbl ), DIMENSION ( : ), ALLOCATABLE :: f,df,fp,p
REAL ( dbl ) :: h, x, res, erf, z, b(3), rmax
Integer :: n, ierr, i, l
TYPE ( radial_grid_type ) :: rg
!deb
!------------------------------------------------------------------------------
CALL timeset ( 'ATOMIC_CODE', 'I', ' ', handle )
IF ( globenv % ionode ) CALL atom_header ( globenv % scr )
!debdebdebdebdeb
n = 100
z=1._dbl
b(1)=1._dbl
b(2)=4._dbl
rmax=10._dbl
CALL init_radial_grid ( rg, "LINEAR", n, z, b, rmax )
! CALL init_radial_grid ( rg, "LOG", n, z, b, rmax )
! CALL init_radial_grid ( rg, "POWER", n, z, b, rmax )
! CALL init_radial_grid ( rg, "RATIONAL", n, z, b, rmax )
! CALL init_radial_grid ( rg, "ALGEBRAIC SQUARE ROOT", n, z, b, rmax )
allocate ( f(0:n+1) )
allocate ( p(0:n+1) )
allocate ( fp(0:n+1) )
allocate ( df(0:n+1) )
do i = 0, n+1
x = rg % r ( i )
f ( i ) = exp ( - x*x )
p ( i ) = x * f ( i )
fp ( i ) = -2._dbl*x*f(i)
enddo
l = 0
call rad_kinetic ( p, p, l, rg, res )
res = res * 4._dbl * pi
write ( 6,*) ' KINETIC ',res
write ( 6,*) ' KINETIC ',1.5_dbl * ( 0.5_dbl*pi )**1.5_dbl
write ( 6,*) ' KINETIC ',res - 1.5_dbl * ( 0.5_dbl*pi )**1.5_dbl
call radial_diff ( n+2, f(0:n+1),df(0:n+1),1,rg,ierr )
call radial_integrate ( n+2,f(0:n+1),rg,res,ierr )
write ( 6,*) ' INTEGRAL ',res
write ( 6,*) ' INTEGRAL ',0.5_dbl * sqrt ( pi ) !* erf (rmax)
write ( 6,*) ' INTEGRAL ',res - 0.5_dbl * sqrt ( pi ) !* erf (rmax)
write ( 6,*) ' DIFF ',sum ( abs ( df - fp ) ) / real ( n )
write ( 6,*) ' DIFF ',minval ( abs ( df - fp ) )
write ( 6,*) ' DIFF ',maxval ( abs ( df - fp ) )
do i = 0, n+1
x = rg % r ( i )
f ( i ) = exp ( - x*x )
fp ( i ) = (4._dbl*x*x-2._dbl)*f(i)
enddo
call radial_diff ( n+2, f(0:n+1),df(0:n+1),2,rg,ierr )
write ( 6,*) ' d2DIFF ',sum ( abs ( df - fp ) ) / real ( n )
write ( 6,*) ' d2DIFF ',minval ( abs ( df - fp ) )
write ( 6,*) ' d2DIFF ',maxval ( abs ( df - fp ) )
do i = 0, n+1
x = rg % r ( i )
f ( i ) = -z* (pi)**(-1.5) * x*x*exp ( - x*x ) * 4._dbl*pi
p ( i ) = -z*erf ( x )
enddo
! call rg_poisson ( f, fp, l, rg, "NUMEROV" )
call rg_poisson ( f, fp, l, rg, "DIFF5P" )
write ( 6,*) ' POISSON ',sum ( abs ( fp - p ) ) / real ( n )
write ( 6,*) ' POISSON ',minval ( abs ( fp - p ) )
write ( 6,*) ' POISSON ',maxval ( abs ( fp - p ) )
do i = 0, n+1
write(99,*) fp(i),p(i)
enddo
deallocate ( f,df,fp,p )
!debdebdebdebdeb
CALL timestop ( zero, handle )
END SUBROUTINE atomic_code
!!*****
!******************************************************************************
END MODULE atom_module
!******************************************************************************

58
src/atoms_input.F
View file

@ -2,6 +2,20 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/atoms_input [1.0] *
!!
!! NAME
!! atoms_input
!!
!! FUNCTION
!!
!! AUTHOR
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE atoms_input
@ -28,21 +42,38 @@ MODULE atoms_input
REAL ( dbl ), POINTER :: v ( :, : )
END TYPE system_type
!!*****
!******************************************************************************
CONTAINS
!!>---------------------------------------------------------------------------!
!! !
!! OPTIONS: either read the section &atoms in the input file or !
!! read the coordinates from the file project_name.dat !
!! SECTION: &atoms ... &end !
!! file: file_name !
!! cell: b11 b12 b13 & !
!! b21 b22 b23 & !
!! b31 b32 b33 !
!! !
!! !
!! !
!!<---------------------------------------------------------------------------!
!******************************************************************************
!!****** atoms_input/read_coord_vel [1.0] *
!!
!! NAME
!! read_coord_vel
!!
!! FUNCTION
!!
!! AUTHOR
!!
!! MODIFICATION HISTORY
!! none
!!
!! OPTIONS
!! either read the section &atoms in the input file or
!! read the coordinates from the file project_name.dat
!!
!! INPUTS
!!I---------------------------------------------------------------------------I
!!I SECTION: &atoms ... &end I
!!I file: file_name I
!!I cell: b11 b12 b13 & I
!!I b21 b22 b23 & I
!!I b31 b32 b33 I
!!I---------------------------------------------------------------------------I
!! SOURCE
!******************************************************************************
SUBROUTINE read_coord_vel ( atype, filen, globenv )
@ -132,6 +163,7 @@ SUBROUTINE read_coord_vel ( atype, filen, globenv )
END SUBROUTINE read_coord_vel
!!*****
!******************************************************************************
SUBROUTINE read_file ( filen, atype )

49
src/band.F
View file

@ -2,6 +2,21 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!------------------------------------------------------------------------------!
!!****** cp2k/band [1.0] *
!!
!! NAME
!! band
!!
!! FUNCTION
!!
!! AUTHOR
!! JGH
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
!
MODULE band
!
@ -28,11 +43,28 @@
REAL (dbl), PARAMETER :: precpar = 1.E-10_dbl, &
delta = 1.E-6_dbl/27.212_dbl
REAL (dbl), PARAMETER :: toll = 1.E-13_dbl
!!*****
!------------------------------------------------------------------------------!
!
CONTAINS
!
!------------------------------------------------------------------------------!
!!****** band/init_band_structure [1.0] *
!!
!! NAME
!! init_band_structure
!!
!! FUNCTION
!!
!! AUTHOR
!! JGH
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
SUBROUTINE init_band_structure(bs,fd,kp)
IMPLICIT NONE
TYPE (band_structure_type), INTENT (INOUT) :: bs
@ -58,7 +90,23 @@
'bs%oc',ns*nk*nspin)
END SUBROUTINE init_band_structure
!!*****
!------------------------------------------------------------------------------!
!!****** band/band_structure_info [1.0] *
!!
!! NAME
!! band_structure_info
!!
!! FUNCTION
!!
!! AUTHOR
!! JGH
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
SUBROUTINE band_structure_info(bs,punit)
IMPLICIT NONE
TYPE (band_structure_type), INTENT (IN) :: bs
@ -127,6 +175,7 @@
WRITE (punit,'(1x,79("-"),/)')
END SUBROUTINE band_structure_info
!!*****
!------------------------------------------------------------------------------!
SUBROUTINE adjust_mu(bs)
!

76
src/brillouin.F
View file

@ -2,6 +2,21 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!------------------------------------------------------------------------------!
!!****** cp2k/brillouin [1.0] *
!!
!! NAME
!! brillouin
!!
!! FUNCTION
!!
!! AUTHOR
!! JGH
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
!
MODULE brillouin
!
@ -29,21 +44,40 @@
REAL (dbl), DIMENSION (:), POINTER :: weight
REAL (dbl), DIMENSION (:,:), POINTER :: xk
END TYPE kpoint_type
!!*****
!------------------------------------------------------------------------------!
!
CONTAINS
!
!!>----------------------------------------------------------------------------!
!! SECTION: &kpoint... &end !
!! !
!! scheme [Gamma, Monkhorst-Pack, MacDonald, General] !
!! { nx ny nz } !
!! { nx ny nz sx sy sz } !
!! { nkpt x1 y1 z1 w1 ... xn yn zn wn } !
!! symmetry [on, off] !
!! wavefunction [real, complex] !
!! !
!!<----------------------------------------------------------------------------!
!------------------------------------------------------------------------------!
!!****** brillouin/kpoint_input [1.0] *
!!
!! NAME
!! kpoint_input
!!
!! FUNCTION
!!
!! AUTHOR
!! JGH
!!
!! MODIFICATION HISTORY
!! none
!!
!! INPUTS
!!I----------------------------------------------------------------------------I
!!I SECTION: &kpoint... &end I
!!I I
!!I scheme [Gamma, Monkhorst-Pack, MacDonald, General] I
!!I { nx ny nz } I
!!I { nx ny nz sx sy sz } I
!!I { nkpt x1 y1 z1 w1 ... xn yn zn wn } I
!!I symmetry [on, off] I
!!I wavefunction [real, complex] I
!!I I
!!I----------------------------------------------------------------------------I
!!
!! SOURCE
!******************************************************************************
SUBROUTINE kpoint_input(kp,inpar)
IMPLICIT NONE
@ -155,7 +189,23 @@
kp%scheme = 'GAMMA'
END IF
END SUBROUTINE kpoint_input
!!*****
!------------------------------------------------------------------------------!
!!****** brillouin/brillouin_info [1.0] *
!!
!! NAME
!! brillouin_info
!!
!! FUNCTION
!!
!! AUTHOR
!! JGH
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
SUBROUTINE brillouin_info(kp,punit)
IMPLICIT NONE
TYPE (kpoint_type), INTENT (IN) :: kp
@ -198,6 +248,10 @@
END DO
END IF
END SUBROUTINE brillouin_info
!!*****
!------------------------------------------------------------------------------!
END MODULE brillouin
!------------------------------------------------------------------------------!

64
src/cntl_input.F
View file

@ -2,6 +2,20 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/cntl_input [1.0] *
!!
!! NAME
!! cntl_input
!!
!! FUNCTION
!!
!! AUTHOR
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE cntl_input
@ -16,22 +30,41 @@ MODULE cntl_input
PRIVATE
PUBLIC :: read_cntl_section
!!*****
!******************************************************************************
CONTAINS
!!>---------------------------------------------------------------------------!
!! SECTION: &cntl ... &end !
!! !
!! simulation [md,debug] !
!! printlevel globenv%print_level !
!! units [kelvin,atomic] !
!! periodic [0,1][0,1][0,1] !
!! Ewald_type [pme_gauss,ewald_gauss] !
!! Ewald_param alpha,[gmax,ns_max,epsilon] !
!! set_file "filename" !
!! input_file "filename" !
!! symmetry [on,off]
!! !
!!<---------------------------------------------------------------------------!
!******************************************************************************
!!****** cntl_input/read_cntl_section [1.0] *
!!
!! NAME
!! read_cntl_section
!!
!! FUNCTION
!!
!! AUTHOR
!!
!! MODIFICATION HISTORY
!! none
!!
!! INPUTS
!!I---------------------------------------------------------------------------I
!!I SECTION: &cntl ... &end I
!!I I
!!I simulation [md,debug] I
!!I printlevel globenv%print_level I
!!I units [kelvin,atomic] I
!!I periodic [0,1][0,1][0,1] I
!!I Ewald_type [pme_gauss,ewald_gauss] I
!!I Ewald_param alpha,[gmax,ns_max,epsilon] I
!!I set_file "filename" I
!!I input_file "filename" I
!!I symmetry [on,off] I
!!I I
!!I---------------------------------------------------------------------------I
!! SOURCE
!******************************************************************************
SUBROUTINE read_cntl_section ( setup, ewald_param, globenv )
@ -173,6 +206,9 @@ SUBROUTINE read_cntl_section ( setup, ewald_param, globenv )
END SUBROUTINE read_cntl_section
!!*****
!******************************************************************************
END MODULE cntl_input
!******************************************************************************

17
src/coefficient_types.F
View file

@ -2,6 +2,20 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****s* cp2k/coefficient_types [1.0] *
!!
!! NAME
!! coefficient_types
!!
!! FUNCTION
!!
!! AUTHOR
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE coefficient_types
@ -42,6 +56,9 @@ MODULE coefficient_types
INTEGER, PARAMETER :: PW_COMPLEXDATA1D = COMPLEXDATA1D
INTEGER, PARAMETER :: PW_REALDATA3D = REALDATA3D
INTEGER, PARAMETER :: PW_COMPLEXDATA3D = COMPLEXDATA3D
!!*****
!******************************************************************************
CONTAINS

17
src/coefficients.F
View file

@ -2,6 +2,20 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/coefficients [1.0] *
!!
!! NAME
!! coefficients
!!
!! FUNCTION
!!
!! AUTHOR
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE coefficients
@ -18,6 +32,9 @@ MODULE coefficients
MODULE PROCEDURE integral_a, integral_ab
END INTERFACE
!!*****
!******************************************************************************
CONTAINS
!******************************************************************************

17
src/constraint.F
View file

@ -2,6 +2,20 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/constraint [1.0] *
!!
!! NAME
!! constraint
!!
!! FUNCTION
!!
!! AUTHOR
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE constraint
@ -23,6 +37,9 @@ MODULE constraint
REAL ( dbl ) :: scale
REAL ( dbl ), DIMENSION (3,3) :: r_rattle
!!*****
!******************************************************************************
CONTAINS

66
src/control_module.F
View file

@ -2,6 +2,21 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/control_module [1.0] *
!!
!! NAME
!! control_module
!!
!! FUNCTION
!!
!! AUTHOR
!! JGH, 14-NOV-2000
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE control_module
@ -20,7 +35,6 @@ MODULE control_module
USE fist_debug, ONLY : fist_debug_control => debug_control
USE tbmd_debug, ONLY : tbmd_debug_control => debug_control
USE tbmd_initialize, ONLY : tbmd_init, tb_get_numel
USE tbmd_global, ONLY : tbatom, tbhop
USE tbmd_input, ONLY : read_tb_hamiltonian, read_tb_hopping_elements
USE cntl_input, ONLY : read_cntl_section
USE force_fields, ONLY : read_force_field_section, ATOMNAMESLENGTH
@ -53,12 +67,31 @@ MODULE control_module
PUBLIC :: control
TYPE ( mdio_parameters_type ) :: mdio
!!*****
!-----------------------------------------------------------------------------!
CONTAINS
!-----------------------------------------------------------------------------!
! CNTL CNTL CNTL CNTL CNTL CNTL CNTL CNTL CNTL CNTL CNTL CNTL CNTL CNTL CNTL !
!-----------------------------------------------------------------------------!
!!****** control_module/control [1.0] *
!!
!! NAME
!! control
!!
!! FUNCTION
!! Controls program flow for FIST and TBMD
!!
!! AUTHOR
!! JGH
!!
!! MODIFICATION HISTORY
!! Harald Forbert (Dec-2000): enable multiple linked lists
!!
!! SOURCE
!******************************************************************************
SUBROUTINE control ( globenv )
@ -97,6 +130,23 @@ SUBROUTINE control ( globenv )
CALL timeset ( 'CONTROL', 'I', ' ', handle1 )
CALL timeset ( 'CNTL_INIT', 'I', ' ', handle2 )
!
struc % ll_data % natom_types = 0
struc % ll_data % list_type = 0
struc % ll_data % counter = 0
struc % ll_data % last_update = 0
struc % ll_data % num_update = 0
struc % ll_data % print_level = 1
NULLIFY ( struc % ll_data % r_last_update )
NULLIFY ( struc % ll_data % rlist_cut )
NULLIFY ( struc % ll_data % rlist_cutsq )
NULLIFY ( struc % ll_data % pp_images )
NULLIFY ( struc % ll_data % pp_ncell )
NULLIFY ( struc % ll_data % pp_startlist_im )
NULLIFY ( struc % ll_data % pp_startlist_nl )
NULLIFY ( struc % ll_data % pp_neighbor )
NULLIFY ( struc % ll_data % pp_startlist_cell )
NULLIFY ( struc % ll_data % pp_cell_ll )
IF ( globenv % program_name == 'FIST' ) THEN
fist = .true.
@ -185,6 +235,7 @@ SUBROUTINE control ( globenv )
! so it also has the right units
CALL get_hinv ( struc % box )
!..initialize molecule_type
CALL initialize_molecule_type ( mol_setup, intra_param, struc % pnode, &
struc % part, struc % molecule, globenv )
@ -210,10 +261,10 @@ SUBROUTINE control ( globenv )
rcut ( :, : ) = inter%potparm ( :, : ) % rcutsq
CALL set_ll_parm ( globenv, simpar % verlet_skin, &
CALL set_ll_parm ( struc%ll_data, globenv, simpar % verlet_skin, &
setup % natom_type, rcut, simpar % n_cell )
CALL set_ll_parm ( globenv, printlevel = globenv % print_level, &
CALL set_ll_parm ( struc%ll_data, globenv, printlevel = globenv % print_level, &
ltype = 'NONBOND' )
DEALLOCATE ( rcut, STAT = isos )
@ -231,7 +282,7 @@ SUBROUTINE control ( globenv )
END IF
!
!..initialize the on-site terms for TB
IF ( tbmd ) CALL tbmd_init ( struc % part, tbatom, tbhop )
IF ( tbmd ) CALL tbmd_init ( struc % part, inter%tbatom, inter%tbhop )
!
!..Symmetry and K-points
!
@ -246,7 +297,7 @@ SUBROUTINE control ( globenv )
!..nota bene: we use the defaults for the elctron temperature and
! spin polarisation
IF ( tbmd ) then
nel = tb_get_numel ( tbatom, charge )
nel = tb_get_numel ( inter%tbatom, charge )
CALL init_fermi_dist(fd,nel)
IF (globenv%ionode) CALL fermi_info(fd,globenv%print_level,globenv%scr)
CALL init_band_structure(bs,fd,kp)
@ -262,10 +313,10 @@ SUBROUTINE control ( globenv )
IF ( fist ) then
CALL fist_debug_control( globenv, ewald_param, struc%part, struc%pnode, &
struc%molecule, struc%box, thermo, inter%potparm, simpar % ensemble )
struc%molecule, struc%box, thermo, inter%potparm, simpar % ensemble, struc%ll_data )
ELSE IF ( tbmd ) THEN
CALL tbmd_debug_control( globenv, ewald_param, struc%part, struc%pnode, &
struc%molecule, struc%box, thermo, inter%potparm, simpar % ensemble )
struc%molecule, struc%box, thermo, inter%potparm, simpar % ensemble , struc%ll_data )
ENDIF
ELSE
@ -303,6 +354,7 @@ SUBROUTINE control ( globenv )
END SUBROUTINE control
!!*****
!-----------------------------------------------------------------------------!
! CNTL CNTL CNTL CNTL CNTL CNTL CNTL CNTL CNTL CNTL CNTL CNTL CNTL CNTL CNTL !
!-----------------------------------------------------------------------------!

25
src/convert_units.F
View file

@ -2,6 +2,20 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/convert_units [1.0] *
!!
!! NAME
!! convert_units
!!
!! FUNCTION
!!
!! AUTHOR
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE convert_units
@ -20,6 +34,9 @@ MODULE convert_units
PRIVATE
PUBLIC :: convert
!!*****
!-----------------------------------------------------------------------------!
CONTAINS
!-----------------------------------------------------------------------------!
@ -145,9 +162,9 @@ SUBROUTINE convert ( units, simpar, part, pstat, box, potparm, &
intra_param % bend_param(:,:,:) % k &
= intra_param % bend_param(:,:,:) % k * units%econv
DO i = 1, 6
intra_param % torsion_param(:,:,:,:) % parm(i) = &
intra_param % torsion_param(:,:,:,:) % parm(i)*units%econv
DO i = 0, SIZE ( intra_param % torsion_param(1,1,1,1) % a ) - 1
intra_param % torsion_param(:,:,:,:) % a(i) = &
intra_param % torsion_param(:,:,:,:) % a(i)*units%econv
END DO
intra_param % constraint_distance ( :, : ) &
= intra_param % constraint_distance ( :, : ) * units%lconv
@ -158,3 +175,5 @@ END SUBROUTINE convert
!******************************************************************************
END MODULE convert_units
!******************************************************************************

75
src/cp2k.F
View file

@ -1,27 +1,44 @@
!-----------------------------------------------------------------------------!
! !
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
! !
! This program is free software; you can redistribute it and/or modify !
! it under the terms of the GNU General Public License as published by !
! the Free Software Foundation; either version 2 of the License, or !
! (at your option) any later version. !
! !
! This program is distributed in the hope that it will be useful, !
! but WITHOUT ANY WARRANTY; without even the implied warranty of !
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the !
! GNU General Public License for more details. !
! !
! You should have received a copy of the GNU General Public License !
! along with this program; if not, write to the Free Software !
! Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. !
! !
!-----------------------------------------------------------------------------!
!
!! Main Program
!! - performs only basic initialisations
!! - calls major modules
!!****** cp2k/cp2k [1.0] *
!!
!! NAME
!! cp2k, The main program
!!
!! COPYRIGHT
!!I---------------------------------------------------------------------------I
!!I I
!!I CP2K: A general program to perform molecular dynamics simulations I
!!I Copyright (C) 2000 CP2K developers group I
!!I I
!!I This program is free software; you can redistribute it and/or modify I
!!I it under the terms of the GNU General Public License as published by I
!!I the Free Software Foundation; either version 2 of the License, or I
!!I (at your option) any later version. I
!!I I
!!I This program is distributed in the hope that it will be useful, I
!!I but WITHOUT ANY WARRANTY; without even the implied warranty of I
!!I MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the I
!!I GNU General Public License for more details. I
!!I I
!!I You should have received a copy of the GNU General Public License I
!!I along with this program; if not, write to the Free Software I
!!I Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. I
!!I I
!!I---------------------------------------------------------------------------I
!!
!! FUNCTION
!! Main Program
!!
!! AUTHOR
!! CJM & JGH
!!
!! MODIFICATION HISTORY
!! none
!!
!! NOTES
!! Main Program
!!
!! SOURCE
!------------------------------------------------------------------------------
PROGRAM cp2k
@ -30,6 +47,7 @@ PROGRAM cp2k
USE control_module, ONLY : control
USE global_types, ONLY : global_environment_type
USE library_tests,ONLY : lib_test
USE atom_module,ONLY : atomic_code
USE kinds, ONLY : dbl, print_kind_info
USE parallel, ONLY : start_parallel, end_parallel
USE physcon, ONLY : print_physcon
@ -56,6 +74,10 @@ PROGRAM cp2k
CALL control ( globenv )
ELSE IF ( globenv % program_name == 'ATOM' ) THEN
CALL atomic_code ( globenv )
ELSE IF ( globenv % program_name == 'TEST' ) THEN
CALL lib_test ( globenv )
@ -66,4 +88,9 @@ PROGRAM cp2k
CALL end_parallel ( )
!------------------------------------------------------------------------------
END PROGRAM cp2k
!------------------------------------------------------------------------------
!!*****

80
src/cp2k_input.F
View file

@ -2,6 +2,22 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/cp2k_input [1.0] *
!!
!! NAME
!! cp2k_input
!!
!! FUNCTION
!! Master input
!!
!! AUTHOR
!! JGH
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE cp2k_input
@ -11,22 +27,45 @@ MODULE cp2k_input
USE parser, ONLY : parser_init, parser_end, read_line, test_next, &
cfield, p_error, get_real, get_int
USE stop_program, ONLY : stop_prg
USE string_utilities, ONLY : uppercase
USE string_utilities, ONLY : uppercase, xstring
USE fft_tools, ONLY : init_fft
PRIVATE
PUBLIC :: read_cp2k_section
!!*****
!******************************************************************************
CONTAINS
!!>---------------------------------------------------------------------------!
!! !
!! SECTION: &cp2k ... &end !
!! !
!! program prg_name !
!! project name !
!! iolevel n !
!! !
!!<---------------------------------------------------------------------------!
!******************************************************************************
!!****** cp2k_input/read_cp2k_section [1.0] *
!!
!! NAME
!! read_cp2k_section
!!
!! FUNCTION
!!
!! AUTHOR
!! JGH
!!
!! MODIFICATION HISTORY
!! 2-Dec-2000 (JGH) added default fft library
!!
!! INPUTS
!!I---------------------------------------------------------------------------I
!!I I
!!I SECTION: &cp2k ... &end I
!!I I
!!I program prg_name I
!!I project name I
!!I iolevel n I
!!I fftlib default_fft_library I
!!I I
!!I---------------------------------------------------------------------------I
!!
!! SOURCE
!******************************************************************************
SUBROUTINE read_cp2k_section ( globenv )
@ -36,7 +75,7 @@ SUBROUTINE read_cp2k_section ( globenv )
TYPE ( global_environment_type ), INTENT ( INOUT ) :: globenv
! Locals
INTEGER :: ierror, ilen, iw, source, allgrp
INTEGER :: ierror, ilen, iw, source, allgrp, i1, i2
CHARACTER ( LEN = 6 ) :: string
CHARACTER ( LEN = 5 ) :: label
@ -44,6 +83,7 @@ SUBROUTINE read_cp2k_section ( globenv )
globenv%program_name = 'NONE'
globenv%project_name = 'project'
globenv%print_level = 0
globenv%default_fft_library = 'FFTSG'
iw = globenv % scr
! ..parse the input section
@ -77,8 +117,14 @@ SUBROUTINE read_cp2k_section ( globenv )
ilen = 0
CALL cfield ( globenv % project_name, ilen )
CASE ( 'IOLEVE')
CASE ( 'IOLEVE' )
globenv%print_level = get_int()
CASE ( 'FFTLIB' )
ilen = 0
CALL cfield ( globenv % default_fft_library, ilen )
CALL uppercase ( globenv % default_fft_library )
END SELECT
CALL read_line()
@ -89,6 +135,11 @@ SUBROUTINE read_cp2k_section ( globenv )
! ..end of parsing the input section
IF ( globenv % program_name == 'QUICKSTEP' ) globenv % program_name = 'QS'
! initialize FFT library
CALL xstring ( globenv % default_fft_library, i1, i2 )
CALL init_fft ( fftlib = globenv % default_fft_library ( i1:i2 ) )
IF ( globenv % ionode ) THEN
! ..write some information to output
@ -98,6 +149,8 @@ SUBROUTINE read_cp2k_section ( globenv )
ADJUSTR ( globenv % input_file_name ( 1:40 ) )
WRITE ( iw, '( A, T41, A )' ) ' CP2K| Project name ', &
ADJUSTR ( globenv % project_name ( 1:40 ) )
WRITE ( iw, '( A, T41, A )' ) ' CP2K| Default FFT library ', &
ADJUSTR ( globenv % default_fft_library ( 1:40 ) )
WRITE ( iw, '( A, T77, I4 )' ) ' CP2K| Global print level ', &
globenv % print_level
WRITE ( iw, '( A, T75, I6 )' ) ' CP2K| Total number of processors ', &
@ -109,6 +162,9 @@ SUBROUTINE read_cp2k_section ( globenv )
END SUBROUTINE read_cp2k_section
!!*****
!******************************************************************************
END MODULE cp2k_input
!******************************************************************************

17
src/dg_rho0s.F
View file

@ -2,6 +2,20 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/dg_rho0s [1.0] *
!!
!! NAME
!! dg_rho0s
!!
!! FUNCTION
!!
!! AUTHOR
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE dg_rho0s
@ -28,6 +42,9 @@ MODULE dg_rho0s
REAL ( dbl ), DIMENSION ( : ), POINTER :: zet
TYPE ( coeff_type ) :: density
END TYPE dg_rho0_type
!!*****
!******************************************************************************
CONTAINS

16
src/dg_types.F
View file

@ -2,6 +2,19 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****s* cp2k/dg_types [1.0] *
!!
!! NAME
!! dg_types
!!
!! FUNCTION
!!
!! AUTHOR
!!
!! MODIFICATION HISTORY
!! none
!!
!******************************************************************************
MODULE dg_types
USE dg_rho0s, ONLY : dg_rho0_type
@ -21,5 +34,8 @@ MODULE dg_types
PUBLIC :: dg_type
!-----------------------------------------------------------------------------!
END MODULE dg_types
!-----------------------------------------------------------------------------!
!!*****

17
src/dgs.F
View file

@ -2,6 +2,20 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/dgs [1.0] *
!!
!! NAME
!! dgs
!!
!! FUNCTION
!!
!! AUTHOR
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE dgs
@ -27,6 +41,9 @@ MODULE dgs
dg_sum_patch_force, dg_get_1x_patch, dg_get_2x_patch, &
dg_get_strucfac, dg_get_rho_tot
!!*****
!******************************************************************************
CONTAINS
!******************************************************************************

34
src/dump.F
View file

@ -2,6 +2,20 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/dump [1.0] *
!!
!! NAME
!! dump
!!
!! FUNCTION
!!
!! AUTHOR
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE dump
@ -14,8 +28,25 @@ MODULE dump
PRIVATE
PUBLIC :: dump_variables
!!*****
!******************************************************************************
CONTAINS
!******************************************************************************
!!****** dump/dump_variables [1.0] *
!!
!! NAME
!! dump_variables
!!
!! FUNCTION
!!
!! AUTHOR
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
SUBROUTINE dump_variables ( struc, dump_file_name, globenv )
@ -57,6 +88,9 @@ SUBROUTINE dump_variables ( struc, dump_file_name, globenv )
END SUBROUTINE dump_variables
!!*****
!******************************************************************************
END MODULE dump
!******************************************************************************

24
src/eigenvalueproblems.F
View file

@ -2,9 +2,24 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!! We are using LAPACK interfaces, so please make sure in IBM/AIX you have
!! the lapack library before essl: "xlf90 ... -llapack -lessl" !!!
!!****** cp2k/eigenvalueproblems [1.0] *
!!
!! NAME
!! eigenvalueproblems
!!
!! FUNCTION
!!
!! AUTHOR
!!
!! MODIFICATION HISTORY
!! none
!!
!! NOTES
!! We are using LAPACK interfaces, so please make sure in IBM/AIX you have
!! the lapack library before essl: "xlf90 ... -llapack -lessl" !!!
!!
!! SOURCE
!******************************************************************************
MODULE eigenvalueproblems
@ -29,6 +44,9 @@ MODULE eigenvalueproblems
MODULE PROCEDURE diagonalise_zhpgv
END INTERFACE
!!*****
!******************************************************************************
CONTAINS
!******************************************************************************

81
src/environment.F
View file

@ -2,6 +2,25 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/environment [1.0] *
!!
!! NAME
!! environment
!!
!! FUNCTION
!! Sets up and terminates the global environment vaiables
!!
!! AUTHOR
!! JGH
!!
!! MODIFICATION HISTORY
!! none
!!
!! USED BY
!! cp2k
!!
!! SOURCE
!******************************************************************************
MODULE environment
@ -15,7 +34,6 @@ MODULE environment
USE timesl, ONLY : walltime, cputime, datum
USE timings, ONLY : timeprint, timeset, timestop
USE util, ONLY : ran2
USE fft_tools, ONLY : init_fft
IMPLICIT NONE
@ -28,9 +46,28 @@ MODULE environment
INTEGER :: my_pid = 0 !! process ID
INTEGER :: my_uid = 0 !! user ID
INTEGER :: handle
!!*****
!******************************************************************************
CONTAINS
!******************************************************************************
!!****** environment/initialisation [1.0] *
!!
!! NAME
!! initialisation
!!
!! FUNCTION
!! Initializes global environment vaiables
!!
!! AUTHOR
!! JGH
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
SUBROUTINE initialisation ( globenv )
@ -42,7 +79,7 @@ SUBROUTINE initialisation ( globenv )
! Locals
CHARACTER ( LEN = 26 ) :: datx
INTEGER :: iw, l, initial_random
INTEGER :: iw, l, initial_random, i1, i2
!------------------------------------------------------------------------------
@ -88,11 +125,25 @@ SUBROUTINE initialisation ( globenv )
! initialize physical constants
CALL init_physcon()
! initialize FFT library
CALL init_fft ( )
END SUBROUTINE initialisation
!!*****
!******************************************************************************
!!****** environment/trailer [1.0] *
!!
!! NAME
!! trailer
!!
!! FUNCTION
!! Writes final timings and banner
!!
!! AUTHOR
!! JGH
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
SUBROUTINE trailer ( globenv )
@ -139,6 +190,23 @@ SUBROUTINE trailer ( globenv )
END SUBROUTINE trailer
!!*****
!******************************************************************************
!!****** environment/set_environment [1.0] *
!!
!! NAME
!! set_environment
!!
!! FUNCTION
!! Calls for UNIX environment variabls
!!
!! AUTHOR
!! JGH
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
SUBROUTINE set_environment ( globenv )
@ -175,6 +243,9 @@ SUBROUTINE set_environment ( globenv )
END SUBROUTINE set_environment
!!*****
!******************************************************************************
END MODULE environment
!******************************************************************************

319
src/erf.F
View file

@ -1,319 +0,0 @@
SUBROUTINE CALERF(ARG,RESULT,JINT)
!------------------------------------------------------------------
!
! This packet evaluates erf(x), erfc(x), and exp(x*x)*erfc(x)
! for a real argument x. It contains three FUNCTION type
! subprograms: ERF, ERFC, and ERFCX (or DERF, DERFC, and DERFCX),
! and one SUBROUTINE type subprogram, CALERF. The calling
! statements for the primary entries are:
!
! Y=ERF(X) (or Y=DERF(X)),
!
! Y=ERFC(X) (or Y=DERFC(X)),
! and
! Y=ERFCX(X) (or Y=DERFCX(X)).
!
! The routine CALERF is intended for internal packet use only,
! all computations within the packet being concentrated in this
! routine. The function subprograms invoke CALERF with the
! statement
!
! CALL CALERF(ARG,RESULT,JINT)
!
! where the parameter usage is as follows
!
! Function Parameters for CALERF
! call ARG Result JINT
!
! ERF(ARG) ANY REAL ARGUMENT ERF(ARG) 0
! ERFC(ARG) ABS(ARG) .LT. XBIG ERFC(ARG) 1
! ERFCX(ARG) XNEG .LT. ARG .LT. XMAX ERFCX(ARG) 2
!
! The main computation evaluates near-minimax approximations
! from "Rational Chebyshev approximations for the error function"
! by W. J. Cody, Math. Comp., 1969, PP. 631-638. This
! transportable program uses rational functions that theoretically
! approximate erf(x) and erfc(x) to at least 18 significant
! decimal digits. The accuracy achieved depends on the arithmetic
! system, the compiler, the intrinsic functions, and proper
! selection of the machine-dependent constants.
!
!*******************************************************************
!*******************************************************************
!
! Explanation of machine-dependent constants
!
! XMIN = the smallest positive floating-point number.
! XINF = the largest positive finite floating-point number.
! XNEG = the largest negative argument acceptable to ERFCX;
! the negative of the solution to the equation
! 2*exp(x*x) = XINF.
! XSMALL = argument below which erf(x) may be represented by
! 2*x/sqrt(pi) and above which x*x will not underflow.
! A conservative value is the largest machine number X
! such that 1.0 + X = 1.0 to machine precision.
! XBIG = largest argument acceptable to ERFC; solution to
! the equation: W(x) * (1-0.5/x**2) = XMIN, where
! W(x) = exp(-x*x)/[x*sqrt(pi)].
! XHUGE = argument above which 1.0 - 1/(2*x*x) = 1.0 to
! machine precision. A conservative value is
! 1/[2*sqrt(XSMALL)]
! XMAX = largest acceptable argument to ERFCX; the minimum
! of XINF and 1/[sqrt(pi)*XMIN].
!
! Approximate values for some important machines are:
!
! XMIN XINF XNEG XSMALL
!
! CDC 7600 (S.P.) 3.13E-294 1.26E+322 -27.220 7.11E-15
! CRAY-1 (S.P.) 4.58E-2467 5.45E+2465 -75.345 7.11E-15
! IEEE (IBM/XT,
! SUN, etc.) (S.P.) 1.18E-38 3.40E+38 -9.382 5.96E-8
! IEEE (IBM/XT,
! SUN, etc.) (D.P.) 2.23D-308 1.79D+308 -26.628 1.11D-16
! IBM 195 (D.P.) 5.40D-79 7.23E+75 -13.190 1.39D-17
! UNIVAC 1108 (D.P.) 2.78D-309 8.98D+307 -26.615 1.73D-18
! VAX D-Format (D.P.) 2.94D-39 1.70D+38 -9.345 1.39D-17
! VAX G-Format (D.P.) 5.56D-309 8.98D+307 -26.615 1.11D-16
!
!
! XBIG XHUGE XMAX
!
! CDC 7600 (S.P.) 25.922 8.39E+6 1.80X+293
! CRAY-1 (S.P.) 75.326 8.39E+6 5.45E+2465
! IEEE (IBM/XT,
! SUN, etc.) (S.P.) 9.194 2.90E+3 4.79E+37
! IEEE (IBM/XT,
! SUN, etc.) (D.P.) 26.543 6.71D+7 2.53D+307
! IBM 195 (D.P.) 13.306 1.90D+8 7.23E+75
! UNIVAC 1108 (D.P.) 26.582 5.37D+8 8.98D+307
! VAX D-Format (D.P.) 9.269 1.90D+8 1.70D+38
! VAX G-Format (D.P.) 26.569 6.71D+7 8.98D+307
!
!*******************************************************************
!*******************************************************************
!
! Error returns
!
! The program returns ERFC = 0 for ARG .GE. XBIG;
!
! ERFCX = XINF for ARG .LT. XNEG;
! and
! ERFCX = 0 for ARG .GE. XMAX.
!
!
! Intrinsic functions required are:
!
! ABS, AINT, EXP
!
!
! Author: W. J. Cody
! Mathematics and Computer Science Division
! Argonne National Laboratory
! Argonne, IL 60439
!
! Latest modification: March 19, 1990
!
!------------------------------------------------------------------
USE KINDS, ONLY : dbl
IMLICIT NONE
INTEGER I,JINT
REAL(dbl), DIMENSION(5) :: A
REAL(dbl), DIMENSION(4) :: B
REAL(dbl), DIMENSION(9) :: C
REAL(dbl), DIMENSION(8) :: D
REAL(dbl), DIMENSION(6) :: P
REAL(dbl), DIMENSION(5) :: Q
REAL(dbl) :: ARG,DEL,FOUR,HALF,ONE,RESULT,SIXTEN,SQRPI, &
TWO,THRESH,X,XBIG,XDEN,XHUGE,XINF,XMAX,XNEG,XNUM,XSMALL, &
Y,YSQ,ZERO
!------------------------------------------------------------------
! Mathematical constants
!------------------------------------------------------------------
REAL(dbl), PARAMETER :: FOUR=4._dbl
REAL(dbl), PARAMETER :: ONE=4._dbl
REAL(dbl), PARAMETER :: HALF=4._dbl
REAL(dbl), PARAMETER :: TWO=4._dbl
REAL(dbl), PARAMETER :: ZERO=4._dbl
REAL(dbl), PARAMETER :: SQRPI=5.6418958354775628695E-1_dbl
REAL(dbl), PARAMETER :: THRESH=0.46875_dbl
REAL(dbl), PARAMETER :: SIXTEN=16._dbl
!------------------------------------------------------------------
! Machine-dependent constants
!------------------------------------------------------------------
REAL(dbl), PARAMETER :: XINF=1.79E308_dbl
REAL(dbl), PARAMETER :: XNEG=-26.628_dbl
REAL(dbl), PARAMETER :: XSMALL=1.11E-16_dbl
REAL(dbl), PARAMETER :: XBIG=25.543_dbl
REAL(dbl), PARAMETER :: XHUGE=6.71E7_dbl
REAL(dbl), PARAMETER :: XMAX=2.53E307_dbl
!------------------------------------------------------------------
! Coefficients for approximation to erf in first interval
!------------------------------------------------------------------
DATA A/3.16112374387056560D00,1.13864154151050156D02,
1 3.77485237685302021D02,3.20937758913846947D03,
2 1.85777706184603153D-1/
DATA B/2.36012909523441209D01,2.44024637934444173D02,
1 1.28261652607737228D03,2.84423683343917062D03/
!------------------------------------------------------------------
! Coefficients for approximation to erfc in second interval
!------------------------------------------------------------------
DATA C/5.64188496988670089D-1,8.88314979438837594D0,
1 6.61191906371416295D01,2.98635138197400131D02,
2 8.81952221241769090D02,1.71204761263407058D03,
3 2.05107837782607147D03,1.23033935479799725D03,
4 2.15311535474403846D-8/
DATA D/1.57449261107098347D01,1.17693950891312499D02,
1 5.37181101862009858D02,1.62138957456669019D03,
2 3.29079923573345963D03,4.36261909014324716D03,
3 3.43936767414372164D03,1.23033935480374942D03/
!------------------------------------------------------------------
! Coefficients for approximation to erfc in third interval
!------------------------------------------------------------------
DATA P/3.05326634961232344D-1,3.60344899949804439D-1,
1 1.25781726111229246D-1,1.60837851487422766D-2,
2 6.58749161529837803D-4,1.63153871373020978D-2/
DATA Q/2.56852019228982242D00,1.87295284992346047D00,
1 5.27905102951428412D-1,6.05183413124413191D-2,
2 2.33520497626869185D-3/
!------------------------------------------------------------------
X = ARG
Y = ABS(X)
IF (Y <= THRESH) THEN
!------------------------------------------------------------------
! Evaluate erf for |X| <= 0.46875
!------------------------------------------------------------------
YSQ = ZERO
IF (Y > XSMALL) YSQ = Y * Y
XNUM = A(5)*YSQ
XDEN = YSQ
DO I = 1, 3
XNUM = (XNUM + A(I)) * YSQ
XDEN = (XDEN + B(I)) * YSQ
END DO
RESULT = X * (XNUM + A(4)) / (XDEN + B(4))
IF (JINT .NE. 0) RESULT = ONE - RESULT
IF (JINT .EQ. 2) RESULT = EXP(YSQ) * RESULT
GO TO 800
!------------------------------------------------------------------
! Evaluate erfc for 0.46875 <= |X| <= 4.0
!------------------------------------------------------------------
ELSE IF (Y <= FOUR) THEN
XNUM = C(9)*Y
XDEN = Y
DO I = 1, 7
XNUM = (XNUM + C(I)) * Y
XDEN = (XDEN + D(I)) * Y
END DO
RESULT = (XNUM + C(8)) / (XDEN + D(8))
IF (JINT /= 2) THEN
YSQ = AINT(Y*SIXTEN)/SIXTEN
DEL = (Y-YSQ)*(Y+YSQ)
RESULT = EXP(-YSQ*YSQ) * EXP(-DEL) * RESULT
END IF
!------------------------------------------------------------------
! Evaluate erfc for |X| > 4.0
!------------------------------------------------------------------
ELSE
RESULT = ZERO
IF (Y >= XBIG) THEN
IF ((JINT /= 2) .OR. (Y >= XMAX)) GO TO 300
IF (Y >= XHUGE) THEN
RESULT = SQRPI / Y
GO TO 300
END IF
END IF
YSQ = ONE / (Y * Y)
XNUM = P(6)*YSQ
XDEN = YSQ
DO I = 1, 4
XNUM = (XNUM + P(I)) * YSQ
XDEN = (XDEN + Q(I)) * YSQ
END DO
RESULT = YSQ *(XNUM + P(5)) / (XDEN + Q(5))
RESULT = (SQRPI - RESULT) / Y
IF (JINT /= 2) THEN
YSQ = AINT(Y*SIXTEN)/SIXTEN
DEL = (Y-YSQ)*(Y+YSQ)
RESULT = EXP(-YSQ*YSQ) * EXP(-DEL) * RESULT
END IF
END IF
!------------------------------------------------------------------
! Fix up for negative argument, erf, etc.
!------------------------------------------------------------------
300 IF (JINT == 0) THEN
RESULT = (HALF - RESULT) + HALF
IF (X < ZERO) RESULT = -RESULT
ELSE IF (JINT .EQ. 1) THEN
IF (X < ZERO) RESULT = TWO - RESULT
ELSE
IF (X < ZERO) THEN
IF (X < XNEG) THEN
RESULT = XINF
ELSE
YSQ = AINT(X*SIXTEN)/SIXTEN
DEL = (X-YSQ)*(X+YSQ)
Y = EXP(YSQ*YSQ) * EXP(DEL)
RESULT = (Y+Y) - RESULT
END IF
END IF
END IF
800 RETURN
!---------- Last card of CALERF ----------
END
FUNCTION ERF(X)
!--------------------------------------------------------------------
!
! This subprogram computes approximate values for erf(x).
! (see comments heading CALERF).
!
! Author/date: W. J. Cody, January 8, 1985
!
!--------------------------------------------------------------------
INTEGER JINT
REAL(dbl) :: ERF, X, RESULT
!------------------------------------------------------------------
JINT = 0
CALL CALERF(X,RESULT,JINT)
ERF = RESULT
!---------- Last card of ERF ----------
END FUNCTION ERF
FUNCTION ERFC(X)
!--------------------------------------------------------------------
!
! This subprogram computes approximate values for erfc(x).
! (see comments heading CALERF).
!
! Author/date: W. J. Cody, January 8, 1985
!
!--------------------------------------------------------------------
INTEGER JINT
REAL(dbl) :: ERFC, X, RESULT
!------------------------------------------------------------------
JINT = 1
CALL CALERF(X,RESULT,JINT)
ERFC = RESULT
C---------- Last card of DERFC ----------
END FUNCTION ERFC
FUNCTION ERFCX(X)
!------------------------------------------------------------------
!
! This subprogram computes approximate values for exp(x*x) * erfc(x).
! (see comments heading CALERF).
!
! Author/date: W. J. Cody, March 30, 1987
!
!------------------------------------------------------------------
INTEGER JINT
REAL(dbl) :: ERFCX, X, RESULT
!------------------------------------------------------------------
JINT = 2
CALL CALERF(X,RESULT,JINT)
DERFCX = RESULT
!---------- Last card of DERFCX ----------
END FUNCTION ERFCX

15
src/erf_fn.F
View file

@ -1,3 +1,17 @@
!******************************************************************************
!!****f* cp2k/CALERF [1.0] *
!!
!! NAME
!! CALERF
!!
!! FUNCTION
!! Calculates error function and complementary error function
!!
!! MODIFICATION HISTORY
!! adapted to F95
!!
!! SOURCE
!******************************************************************************
SUBROUTINE CALERF(ARG,RESULT,JINT)
@ -346,3 +360,4 @@ FUNCTION ERFCX(X)
RETURN
!---------- Last card of DERFCX ----------
END
!!*****

17
src/ewald_parameters_types.F
View file

@ -2,6 +2,20 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****s* cp2k/ewald_parameters_types [1.0] *
!!
!! NAME
!! ewald_parameters_types
!!
!! FUNCTION
!!
!! AUTHOR
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE ewald_parameters_types
@ -24,3 +38,6 @@ MODULE ewald_parameters_types
PUBLIC :: ewald_parameters_type
END MODULE ewald_parameters_types
!******************************************************************************
!!*****

17
src/ewalds.F
View file

@ -2,6 +2,20 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/ewalds [1.0] *
!!
!! NAME
!! ewalds
!!
!! FUNCTION
!!
!! AUTHOR
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE ewalds
@ -36,6 +50,9 @@ MODULE ewalds
INTEGER :: ewald_grp
!!*****
!******************************************************************************
CONTAINS
!******************************************************************************

50
src/fermi.F
View file

@ -2,6 +2,21 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!------------------------------------------------------------------------------!
!!****** cp2k/fermi [1.0] *
!!
!! NAME
!! fermi
!!
!! FUNCTION
!!
!! AUTHOR
!! JGH
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
!
MODULE fermi
!
@ -25,11 +40,27 @@
INTEGER :: multiplicity
INTEGER :: nstate, na, nb
END TYPE fermi_distribution_type
!!*****
!------------------------------------------------------------------------------!
!
CONTAINS
!
!------------------------------------------------------------------------------!
!!****** fermi/init_fermi_dist [1.0] *
!!
!! NAME
!! init_fermi_dist
!!
!! FUNCTION
!!
!! AUTHOR
!! JGH
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
SUBROUTINE init_fermi_dist(fdist,nel,etemp,lsd,mult,nstate)
IMPLICIT NONE
TYPE (fermi_distribution_type), INTENT (INOUT) :: fdist
@ -130,7 +161,23 @@
END IF
END IF
END SUBROUTINE init_fermi_dist
!!*****
!------------------------------------------------------------------------------!
!!****** fermi/fermi_info [1.0] *
!!
!! NAME
!! fermi_info
!!
!! FUNCTION
!!
!! AUTHOR
!! JGH
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
SUBROUTINE fermi_info(fdist,plevel,punit)
IMPLICIT NONE
TYPE (fermi_distribution_type), INTENT (IN) :: fdist
@ -171,6 +218,9 @@
END IF
END IF
END SUBROUTINE fermi_info
!!*****
!------------------------------------------------------------------------------!
END MODULE fermi
!------------------------------------------------------------------------------!

118
src/fft_tools.F
View file

@ -2,19 +2,35 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!
! How to add a new FFT library:
! - create a new interface library : fftXX_lib with the entries
! fft3d and fft_get_lengths, see fftw_lib for a template
! - add in this file the entries to the new library; in each
! subroutine there will be an additional CASE
!
!!****** cp2k/fft_tools [1.0] *
!!
!! NAME
!! fft_tools
!!
!! FUNCTION
!!
!! AUTHOR
!! JGH
!!
!! MODIFICATION HISTORY
!! none
!!
!! NOTES
!! How to add a new FFT library:
!! - create a new interface library : fftXX_lib with the entries
!! fft3d and fft_get_lengths, see fftw_lib for a template
!! - add in this file the entries to the new library; in each
!! subroutine there will be an additional CASE
!!
!! SOURCE
!******************************************************************************
MODULE fft_tools
USE kinds, ONLY: dbl, sgl
USE stop_program, ONLY : stop_prg
USE fftessl_lib, ONLY : i_fft3d => fft3d, &
i_fft_get_lengths => fft_get_lengths
USE fftw_lib, ONLY : w_fft3d => fft3d, &
w_fft_get_lengths => fft_get_lengths
USE fftsg_lib, ONLY : sg_fft3d => fft3d, &
@ -35,10 +51,31 @@ MODULE fft_tools
INTEGER :: fft_type = 0
!!*****
!******************************************************************************
CONTAINS
!******************************************************************************
!!****** fft_tools/init_fft [1.0] *
!!
!! NAME
!! init_fft
!!
!! FUNCTION
!!
!! AUTHOR
!! JGH
!!
!! MODIFICATION HISTORY
!! 30-Nov-2000 (JGH) ESSL FFT Library added
!!
!! NOTES
!! Label 1 : FFTsg
!! Label 2 : FFTw
!! Label 3 : FFT ESSL
!!
!! SOURCE
!******************************************************************************
SUBROUTINE init_fft ( fftlib, fftnum )
@ -53,11 +90,13 @@ SUBROUTINE init_fft ( fftlib, fftnum )
IF ( PRESENT ( fftlib ) ) THEN
SELECT CASE ( fftlib )
CASE DEFAULT
CALL stop_prg ("init_fft","Unknown FFT library")
CALL stop_prg ("init_fft","Unknown FFT library",fftlib)
CASE ( "FFTSG" )
fft_type = 1
CASE ( "FFTW" )
fft_type = 2
CASE ( "FFTESSL" )
fft_type = 3
END SELECT
ELSE IF ( PRESENT ( fftnum ) ) THEN
fft_type = fftnum
@ -82,6 +121,22 @@ SUBROUTINE init_fft ( fftlib, fftnum )
END SUBROUTINE init_fft
!!*****
!******************************************************************************
!!****** fft_tools/get_fft_library [1.0] *
!!
!! NAME
!! get_fft_library
!!
!! FUNCTION
!!
!! AUTHOR
!! JGH
!!
!! MODIFICATION HISTORY
!! 30-Nov-2000 (JGH) ESSL FFT Library added
!!
!! SOURCE
!******************************************************************************
SUBROUTINE get_fft_library ( fft_handle, library )
@ -99,10 +154,29 @@ SUBROUTINE get_fft_library ( fft_handle, library )
CASE ( 2 )
library = " FFTw library initialized "
fft_handle = 2
CASE ( 3 )
library = " ESSL FFT library initialized "
fft_handle = 3
END SELECT
END SUBROUTINE get_fft_library
!!*****
!******************************************************************************
!!****** fft_tools/fft_radix_operations [1.0] *
!!
!! NAME
!! fft_radix_operations
!!
!! FUNCTION
!!
!! AUTHOR
!! Ari Seitsonen
!!
!! MODIFICATION HISTORY
!! new library structure (JGH)
!!
!! SOURCE
!******************************************************************************
SUBROUTINE fft_radix_operations ( radix_in, radix_out, operation )
@ -125,6 +199,8 @@ SUBROUTINE fft_radix_operations ( radix_in, radix_out, operation )
CALL sg_fft_get_lengths ( data, max_length = 1024 )
CASE ( 2 )
CALL w_fft_get_lengths ( data, max_length = 1024 )
CASE ( 3 )
CALL i_fft_get_lengths ( data, max_length = 1024 )
CASE DEFAULT
CALL stop_prg ("fft3d","Unknown FFT library")
END SELECT
@ -176,6 +252,23 @@ SUBROUTINE fft_radix_operations ( radix_in, radix_out, operation )
END SUBROUTINE fft_radix_operations
!!*****
!******************************************************************************
!!****** fft_tools/fft3d [1.0] *
!!
!! NAME
!! fft3d
!!
!! FUNCTION
!! Calls the 3D-FFT function from the initialized library
!!
!! AUTHOR
!! JGH
!!
!! MODIFICATION HISTORY
!! 30-Nov-2000 (JGH) ESSL FFT Library added
!!
!! SOURCE
!******************************************************************************
SUBROUTINE fft3d ( fsign, n, zg, zg_out, scale, status )
@ -214,6 +307,12 @@ SUBROUTINE fft3d ( fsign, n, zg, zg_out, scale, status )
ELSE
CALL w_fft3d ( sign, norm, n, zg )
END IF
CASE ( 3 )
IF ( PRESENT ( zg_out ) ) THEN
CALL i_fft3d ( sign, norm, n, zg, zg_out )
ELSE
CALL i_fft3d ( sign, norm, n, zg )
END IF
CASE DEFAULT
CALL stop_prg ("fft3d","Unknown FFT library")
END SELECT
@ -228,6 +327,7 @@ SUBROUTINE fft3d ( fsign, n, zg, zg_out, scale, status )
END SUBROUTINE fft3d
!!*****
!******************************************************************************
END MODULE fft_tools

259
src/fftessl_lib.F Normal file
View file

@ -0,0 +1,259 @@
!-----------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/fftessl_lib [1.0] *
!!
!! NAME
!! fftessl_lib
!!
!! FUNCTION
!! Interface to FFT routines in ESSL library
!!
!! AUTHOR
!! JGH 30-Nov-2000
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE fftessl_lib
USE kinds, ONLY: dbl
USE stop_program, ONLY : stop_memory, stop_prg
USE util, ONLY : sort
IMPLICIT NONE
PRIVATE
PUBLIC :: fft3d
PUBLIC :: fft_get_lengths
!!*****
!******************************************************************************
CONTAINS
!******************************************************************************
!!****** fftessl_lib/fft_get_lengths [1.0] *
!!
!! NAME
!! fft_get_lengths
!!
!! FUNCTION
!! Give the allowed lengths of FFT's '''
!!
!! AUTHOR
!! JGH
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
SUBROUTINE fft_get_lengths ( data, max_length )
IMPLICIT NONE
! Arguments
INTEGER, INTENT ( IN ) :: max_length
INTEGER, DIMENSION ( : ), POINTER :: data
! Locals
INTEGER :: iloc
INTEGER, DIMENSION ( : ), ALLOCATABLE :: idx
INTEGER :: h, i, j, k, m, number, ndata, nmax, allocstat, maxn
INTEGER :: maxn_twos, maxn_threes, maxn_fives
INTEGER :: maxn_sevens, maxn_elevens, maxn_thirteens
!------------------------------------------------------------------------------
! compute ndata
!! FFTW can do arbitrary(?) lenghts, maybe you want to limit them to some
!! powers of small prime numbers though...
maxn_twos = 15
maxn_threes = 2
maxn_fives = 1
maxn_sevens = 1
maxn_elevens = 1
maxn_thirteens = 0
maxn = MIN ( max_length, 37748736 )
ndata = 0
DO h = 0, maxn_twos
nmax = HUGE(0) / 2**h
DO i = 0, maxn_threes
DO j = 0, maxn_fives
DO k = 0, maxn_sevens
DO m = 0, maxn_elevens
number = (3**i) * (5**j) * (7**k) * (11**m)
IF ( number > nmax ) CYCLE
number = number * 2 ** h
IF ( number >= maxn ) CYCLE
ndata = ndata + 1
END DO
END DO
END DO
END DO
END DO
ALLOCATE ( data ( ndata ), idx ( ndata ), STAT = allocstat )
IF ( allocstat /= 0 ) THEN
CALL stop_memory ( "fft_get_lengths", "data, idx", 2*ndata )
END IF
ndata = 0
data ( : ) = 0
DO h = 0, maxn_twos
nmax = HUGE(0) / 2**h
DO i = 0, maxn_threes
DO j = 0, maxn_fives
DO k = 0, maxn_sevens
DO m = 0, maxn_elevens
number = (3**i) * (5**j) * (7**k) * (11**m)
IF ( number > nmax ) CYCLE
number = number * 2 ** h
IF ( number >= maxn ) CYCLE
ndata = ndata + 1
data ( ndata ) = number
END DO
END DO
END DO
END DO
END DO
CALL sort ( data, ndata, idx )
DEALLOCATE ( idx, STAT = allocstat )
IF ( allocstat /= 0 ) THEN
CALL stop_memory ( "fft_get_lengths", "idx" )
END IF
END SUBROUTINE fft_get_lengths
!!*****
!******************************************************************************
!!****** fftessl_lib/fft_3d [1.0] *
!!
!! NAME
!! fft_3d
!!
!! FUNCTION
!! Routine with wrapper for all 3dfft call:
!! Does transform with exp(+ig.r*sign):
!!
!! AUTHOR
!! JGH
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
SUBROUTINE fft3d ( fsign, scale, n, zg, zg_out )
IMPLICIT NONE
! Arguments
INTEGER, INTENT ( INOUT ) :: fsign
REAL ( dbl ), INTENT ( IN ), OPTIONAL :: scale
INTEGER, DIMENSION ( : ), INTENT ( IN ) :: n
COMPLEX ( dbl ), DIMENSION(:,:,:), INTENT ( INOUT ) :: zg
COMPLEX ( dbl ), DIMENSION(:,:,:), INTENT ( INOUT ), OPTIONAL :: zg_out
! Locals
INTEGER :: sign_fft, ldx, ldy, ldz, ldox, ldoy, ldoz, ierr
INTEGER :: nx, ny, nz, na1, na2, isos, naux
COMPLEX ( dbl ), DIMENSION(:), ALLOCATABLE :: aux
LOGICAL :: fft_in_place
!------------------------------------------------------------------------------
IF ( PRESENT ( zg_out ) ) THEN
fft_in_place = .false.
ELSE
fft_in_place = .true.
END IF
sign_fft = fsign
nx = n ( 1 )
ny = n ( 2 )
nz = n ( 3 )
ldx = SIZE ( zg (:,1,1) )
ldy = SIZE ( zg (1,:,1) )
ldz = SIZE ( zg (1,1,:) )
IF( MAX ( ny, nz ) < 252 ) THEN
IF( nx <= 2048 ) THEN
naux = 60000
ELSE
naux = 60000 + NINT( 4.56_dbl * nx )
END IF
ELSE
IF( nx <= 2048 ) THEN
na1 = 60000 + ( 2 * ny + 256 ) * ( MIN ( 64, nx ) + 4.56 )
na2 = 60000 + ( 2 * nz + 256 ) * ( MIN ( 64, nx * ny ) + 4.56 )
ELSE
na1 = 60000 + 4.56 * nx + ( 2 * ny + 256 ) * ( MIN ( 64, nx ) + 4.56 )
naux = 60000 + 4.56 * nx + ( 2 * nz + 256 ) * ( MIN ( 64, nx * ny ) + 4.56 )
END IF
IF( ny >= 252 .AND. nz < 252 ) THEN
naux=na1
ELSE IF( ny < 252 .AND. nz >= 252 ) THEN
naux = na2
ELSE
naux = MAX ( na1, na2 )
END IF
END IF
#if defined ( __FFTESSL )
ALLOCATE ( aux ( naux ), STAT = isos )
IF ( isos /= 0 ) CALL stop_memory ( 'fft3d','aux',naux )
! sign_fft = fsign
IF ( fft_in_place ) THEN
CALL dcft3(zg,ldx,ldx*ldy,zg,ldx,ldx*ldy,nx,ny,nz,sign_fft,scale,aux,naux)
ELSE
ldox = SIZE ( zg_out (:,1,1) )
ldoy = SIZE ( zg_out (1,:,1) )
ldoz = SIZE ( zg_out (1,1,:) )
CALL dcft3(zg,ldx,ldx*ldy,zg_out,ldox,ldox*ldoz,nx,ny,nz,sign_fft,scale,aux,naux)
END IF
DEALLOCATE ( aux, STAT = isos )
IF ( isos /= 0 ) CALL stop_memory ( 'fft3d','aux' )
#else
fsign = 0
#endif
END SUBROUTINE fft3d
!!*****
!******************************************************************************
END MODULE fftessl_lib
!******************************************************************************

56
src/fftsg_lib.F
View file

@ -2,6 +2,21 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/fftsg_lib [1.0] *
!!
!! NAME
!! fftsg_lib
!!
!! FUNCTION
!!
!! AUTHOR
!! JGH
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE fftsg_lib
@ -14,13 +29,28 @@ MODULE fftsg_lib
PUBLIC :: fft3d
PUBLIC :: fft_get_lengths
!!*****
!******************************************************************************
CONTAINS
!******************************************************************************
!! Give the allowed lengths of FFT's '''
!!****** fftsg_lib/fft_get_lengths [1.0] *
!!
!! NAME
!! fft_get_lengths
!!
!! FUNCTION
!! Give the allowed lengths of FFT's '''
!!
!! AUTHOR
!! JGH
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
SUBROUTINE fft_get_lengths ( data, max_length )
@ -58,10 +88,25 @@ SUBROUTINE fft_get_lengths ( data, max_length )
END SUBROUTINE fft_get_lengths
!!*****
!******************************************************************************
!!****** fftsg_lib/fft_3d [1.0] *
!!
!! NAME
!! fft_3d
!!
!! FUNCTION
!! Routine with wrapper for all 3dfft call:
!! Does transform with exp(+ig.r*sign):
!!
!! AUTHOR
!! JGH
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
! routine with wrapper for all fft calls:
! Does transform with exp(+ig.r*sign):
SUBROUTINE fft3d ( fsign, scale, n, zg, zg_out )
@ -148,6 +193,7 @@ SUBROUTINE fft3d ( fsign, scale, n, zg, zg_out )
END SUBROUTINE fft3d
!!*****
!******************************************************************************
END MODULE fftsg_lib

56
src/fftw_lib.F
View file

@ -2,6 +2,21 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/fftw_lib [1.0] *
!!
!! NAME
!! fftw_lib
!!
!! FUNCTION
!!
!! AUTHOR
!! JGH
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE fftw_lib
@ -45,13 +60,28 @@ MODULE fftw_lib
parameter(FFTW_SCRAMBLED_INPUT=8192)
parameter(FFTW_SCRAMBLED_OUTPUT=16384)
!!*****
!******************************************************************************
CONTAINS
!******************************************************************************
!! Give the allowed lengths of FFT's '''
!!****** fftw_lib/fft_get_lengths [1.0] *
!!
!! NAME
!! fft_get_lengths
!!
!! FUNCTION
!! Give the allowed lengths of FFT's '''
!!
!! AUTHOR
!! JGH
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
SUBROUTINE fft_get_lengths ( data, max_length )
@ -140,10 +170,25 @@ SUBROUTINE fft_get_lengths ( data, max_length )
END SUBROUTINE fft_get_lengths
!!*****
!******************************************************************************
!!****** fftw_lib/fft3d [1.0] *
!!
!! NAME
!! fft3d
!!
!! FUNCTION
!! Routine with wrapper for 3dfft calls:
!! Does transform with exp(+ig.r*sign):
!!
!! AUTHOR
!! JGH
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
! routine with wrapper for all fft calls:
! Does transform with exp(+ig.r*sign):
SUBROUTINE fft3d ( fsign, scale, n, zg, zg_out )
@ -306,6 +351,7 @@ SUBROUTINE fft3d ( fsign, scale, n, zg, zg_out )
END SUBROUTINE fft3d
!!*****
!******************************************************************************
END MODULE fftw_lib

330
src/fist.F
View file

@ -1,330 +0,0 @@
!-----------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
MODULE fist
USE atoms_input, ONLY : read_coord_vel, system_type
USE convert_units, ONLY : convert
USE dump, ONLY : dump_variables
USE ewalds, ONLY : ewald_print, ewald_correction
USE ewald_parameters_types, ONLY : ewald_parameters_type
USE header, ONLY : fist_header
USE kinds, ONLY : dbl
USE fist_debug, ONLY : debug_control
USE fist_input, ONLY : read_fist_section
USE fist_global, ONLY : fistpar
USE force_fields, ONLY : read_force_field_section, ATOMNAMESLENGTH
USE initialize_extended_types, ONLY : initialize_extended_type
USE initialize_molecule_types, ONLY : initialize_molecule_type
USE initialize_particle_types, ONLY : initialize_particle_type
USE integrator, ONLY : velocity_verlet, force, set_energy_parm, energy, &
set_integrator
USE input_types, ONLY : setup_parameters_type
USE linklist_control, ONLY : set_ll_parm
USE mathconstants, ONLY : zero
USE md, ONLY : read_md_section, simulation_parameters_type, &
initialize_velocities, thermodynamic_type, mdio_parameters_type
USE molecule_input, ONLY : read_molecule_section, read_setup_section
USE molecule_types, ONLY : molecule_type, intra_parameters_type
USE nose, ONLY : extended_parameters_type
USE pair_potential, ONLY : potentialparm_type, &
potential_f, spline_nonbond_control
USE particle_types, ONLY : particle_prop_type, particle_type
USE simulation_cell, ONLY : cell_type, get_hinv
USE stop_program, ONLY : stop_prg, stop_memory
USE structure_types, ONLY : structure_type
USE timings, ONLY : timeset, timestop
USE unit, ONLY : unit_convert_type, set_units
USE util, ONLY : close_unit, get_share
IMPLICIT NONE
PRIVATE
PUBLIC :: fist_main
TYPE ( mdio_parameters_type ) :: mdio
CONTAINS
!-----------------------------------------------------------------------------!
! FIST FIST FIST FIST FIST FIST FIST FIST FIST FIST FIST FIST FIST FIST FIST !
!-----------------------------------------------------------------------------!
SUBROUTINE fist_main
IMPLICIT NONE
! Locals
REAL ( dbl ) :: cons, ecut, qi, qj
REAL ( dbl ), DIMENSION ( :, : ), ALLOCATABLE :: rcut
INTEGER :: handle, itimes, isos, para_comm_cart
INTEGER :: iat, jat
CHARACTER ( LEN = ATOMNAMESLENGTH ), DIMENSION ( : ), POINTER :: atom_names
CHARACTER ( LEN = 20 ) :: set_fn
TYPE ( potentialparm_type ), DIMENSION ( :, : ), POINTER :: potparm
TYPE ( particle_prop_type ), DIMENSION ( : ), POINTER :: pstat
TYPE ( molecule_type ), DIMENSION ( : ), POINTER :: mol_setup
TYPE ( unit_convert_type ) :: units
TYPE ( simulation_parameters_type ) :: simpar
TYPE ( structure_type ) :: struc
TYPE ( extended_parameters_type ) :: nhcp
TYPE ( thermodynamic_type ) :: thermo
TYPE ( system_type ) :: ainp
TYPE ( ewald_parameters_type ) :: ewald_param
TYPE ( setup_parameters_type ) :: setup
TYPE ( intra_parameters_type ) :: intra_param
!------------------------------------------------------------------------------
CALL timeset ( 'FIST', 'I', ' ', handle )
IF ( fistpar % ionode ) CALL fist_header ( fistpar % scr )
!..read fist section
CALL read_fist_section ( setup, ewald_param, fistpar )
! read from the setup and molecule section of *.set
set_fn = setup % set_file_name
CALL read_setup_section ( mol_setup, set_fn, fistpar )
CALL read_molecule_section ( mol_setup, set_fn, fistpar )
! read force_field information for classical MD
CALL read_force_field_section ( setup, mol_setup, set_fn, &
intra_param, potparm, atom_names, pstat, fistpar )
!..read the input of the molecular dynamics section
CALL read_md_section ( simpar, fistpar, mdio )
simpar % program = 'FIST'
!..read atomic coordinates, velocities (optional) and the simulation box
ainp % rtype = simpar % read_type
CALL read_coord_vel ( ainp, setup % input_file_name, fistpar )
!..initialize box, perd
struc % box % hmat = ainp % box
struc % box % perd = setup % perd
!..initialize working units
CALL set_units ( setup % unit_type, units )
CALL set_energy_parm ( units % pconv, units % econv, units % l_label, &
units % vol_label, units % e_label, units % pres_label, &
units % temp_label, units % angl_label )
!..allocate memory for atoms and molecules
CALL allocmem ( ainp, mol_setup, struc )
!..initialize particle_type
CALL initialize_particle_type ( atom_names, simpar, mol_setup, &
ainp, pstat, struc % part )
!..convert the units
CALL convert ( units = units, simpar = simpar, &
part = struc % part, pstat = pstat, box = struc % box, &
potparm = potparm, intra_param = intra_param, &
ewald_param = ewald_param )
!..calculate the inverse box matrix now after the unit conversion,
! so it also has the right units
CALL get_hinv ( struc % box )
!..initialize molecule_type
CALL initialize_molecule_type ( mol_setup, intra_param, struc % pnode, &
struc % part, struc % molecule, fistpar )
!..initialize extended_parameters_type and get number of degrees of freedom
CALL initialize_extended_type ( struc % box, simpar, &
struc % molecule, mol_setup, nhcp, fistpar )
! initialize velocities if needed
IF ( simpar % read_type == 'POS' ) THEN
CALL initialize_velocities ( simpar, struc % part, fistpar )
END IF
!...initialize splines
potparm ( :, : ) % energy_cutoff = 0.0_dbl
potparm ( :, : ) % e_cutoff_coul = 0.0_dbl
CALL spline_nonbond_control ( potparm, pstat, 5000, ewald_param )
!..set linklist control parameters
ALLOCATE ( rcut ( setup % natom_type, setup % natom_type ), STAT = isos )
IF ( isos /=0 ) CALL stop_memory ( 'fist', 'rcut', 0 )
rcut ( :, : ) = potparm ( :, : ) % rcutsq
CALL set_ll_parm ( fistpar, simpar % verlet_skin, &
setup % natom_type, rcut, simpar % n_cell )
CALL set_ll_parm ( fistpar, printlevel = fistpar % print_level, &
ltype = 'NONBOND' )
DEALLOCATE ( rcut, STAT = isos )
IF ( isos /= 0 ) CALL stop_memory ( 'fist', 'rcut' )
!..deallocate arrays needed for atom input
IF ( ASSOCIATED ( ainp % c ) ) THEN
DEALLOCATE ( ainp % c, STAT = isos )
IF ( isos /= 0 ) CALL stop_memory ( 'fist', 'ainp%c' )
END IF
IF ( ASSOCIATED ( ainp % v ) ) THEN
DEALLOCATE ( ainp % v, STAT = isos )
IF ( isos /= 0 ) CALL stop_memory ( 'fist', 'ainp%v' )
END IF
!..initialize integrator
CALL set_integrator (fistpar % ionode, fistpar % group, fistpar % scr, mdio )
IF ( setup % run_type == 'DEBUG' ) THEN
!..debug the forces
CALL debug_control( fistpar, ewald_param, struc%part, struc%pnode, &
struc%molecule, struc%box, thermo, potparm, simpar % ensemble )
ELSE
!..MD
itimes = 0
CALL force ( struc, potparm, thermo, simpar, ewald_param )
IF ( fistpar % ionode .AND. ewald_param % ewald_type /= 'NONE' ) &
CALL ewald_print ( fistpar % scr, thermo, struc % box, &
units % e_label )
CALL energy ( itimes, cons, simpar, struc, thermo, nhcp )
DO itimes = 1, simpar % nsteps
CALL velocity_verlet ( itimes, cons, simpar, potparm, thermo, &
struc, ewald_param, nhcp )
IF ( MOD ( itimes, mdio % idump ) == 0 ) &
CALL dump_variables ( struc, mdio % dump_file_name, fistpar )
END DO
CALL dump_variables ( struc, mdio % dump_file_name, fistpar )
IF ( fistpar % ionode ) CALL close_unit ( 10, 99 )
END IF
!..deallocate memory for atoms and molecules
CALL deallocmem ( struc )
CALL timestop ( zero, handle )
END SUBROUTINE fist_main
!-----------------------------------------------------------------------------!
! FIST FIST FIST FIST FIST FIST FIST FIST FIST FIST FIST FIST FIST FIST FIST !
!-----------------------------------------------------------------------------!
SUBROUTINE allocmem ( ainp, mol_setup, struc )
IMPLICIT NONE
! Arguments
TYPE ( molecule_type ), DIMENSION ( : ), INTENT ( IN ) :: mol_setup
TYPE ( structure_type ), INTENT ( INOUT ) :: struc
TYPE ( system_type ), INTENT ( IN ) :: ainp
! Locals
INTEGER :: iw, natoms, nnodes, nmol, nmoltype, ios, iat, i
!------------------------------------------------------------------------------
struc % name = 'FIST MOLECULAR SYSTEM'
IF ( fistpar % num_pe == 1 ) THEN
natoms = SIZE ( ainp % c ( 1, : ) )
ALLOCATE ( struc % part ( 1:natoms ), STAT = ios )
IF ( ios /= 0 ) CALL stop_memory ( 'fist', 'part', natoms )
ALLOCATE ( struc % pnode ( 1:natoms ), STAT = ios )
IF ( ios /= 0 ) CALL stop_memory ( 'fist', 'pnode', natoms )
nmol = SUM ( mol_setup ( : ) % num_mol )
ALLOCATE ( struc % molecule ( 1:nmol ), STAT = ios )
IF ( ios /= 0 ) CALL stop_memory ( 'fist', 'molecule', nmol )
IF ( fistpar % ionode .AND. fistpar % print_level > 3 ) THEN
iw = fistpar % scr
WRITE ( iw, '( A )' )
WRITE ( iw, '( A, T71, I10 )' ) &
' FIST| Number of allocated particles ', natoms
WRITE ( iw,'( A, T71, I10 )' ) &
' FIST| Number of allocated particle nodes ', natoms
WRITE ( iw, '( A, T71, I10 )' ) &
' FIST| Number of allocated molecules ', nmol
WRITE ( iw, '( A )' )
END IF
ELSE
!..replicated data
natoms = SIZE ( ainp % c ( 1, : ) )
ALLOCATE ( struc % part ( 1:natoms ), STAT = ios )
IF ( ios /= 0 ) CALL stop_memory ( 'fist', 'part', natoms )
nmoltype = SIZE ( mol_setup )
nmol = 0
nnodes = 0
DO i = 1, nmoltype
nmol = nmol + get_share ( mol_setup ( i ) % num_mol, &
fistpar % num_pe, fistpar % mepos )
nnodes = nnodes + nmol * mol_setup ( i ) % molpar % natom
END DO
ALLOCATE ( struc % molecule ( 1:nmol ), STAT = ios )
IF ( ios /= 0 ) CALL stop_memory ( 'fist', 'molecule' , nmol )
ALLOCATE ( struc % pnode ( 1:nnodes ), STAT = ios )
IF ( ios /= 0 ) CALL stop_memory ( 'fist', 'pnode', nnodes )
IF ( fistpar % ionode .AND. fistpar % print_level > 3 ) THEN
iw = fistpar % scr
WRITE ( iw, '( A )' )
WRITE ( iw, '( A, T71, I10 )' ) &
' FIST| Number of allocated particles ', natoms
WRITE ( iw, '( A, T71, I10 )' ) &
' FIST| Number of allocated particle nodes ', nnodes
WRITE ( iw, '( A, I5, T71, I10 )' ) &
' FIST| Number of allocated molecules on processor ', &
fistpar % mepos, nmol
WRITE ( iw, '( A )' )
END IF
END IF
END SUBROUTINE allocmem
!******************************************************************************
SUBROUTINE deallocmem ( struc )
IMPLICIT NONE
! Arguments
TYPE ( structure_type ), INTENT ( INOUT ) :: struc
! Locals
INTEGER :: ios
!------------------------------------------------------------------------------
DEALLOCATE ( struc % part, STAT = ios )
IF ( ios /= 0 ) CALL stop_memory ( 'fist', 'part' )
DEALLOCATE ( struc % pnode, STAT = ios )
IF ( ios /= 0 ) CALL stop_memory ( 'fist', 'pnode' )
DEALLOCATE ( struc % molecule, STAT = ios )
IF ( ios /= 0 ) CALL stop_memory ( 'fist', 'molecule' )
END SUBROUTINE deallocmem
!******************************************************************************
END MODULE fist

132
src/fist_debug.F
View file

@ -2,6 +2,21 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/fist_debug [1.0] *
!!
!! NAME
!! fist_debug
!!
!! FUNCTION
!!
!! AUTHOR
!! CJM
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE fist_debug
@ -23,20 +38,39 @@ MODULE fist_debug
USE pw_grids, ONLY : pw_find_cutoff, pw_grid_setup
USE simulation_cell, ONLY : cell_type
USE stop_program, ONLY : stop_memory, stop_prg
USE linklist_types, ONLY : linklist_internal_data_type
IMPLICIT NONE
PRIVATE
PUBLIC :: debug_control
!!*****
!******************************************************************************
CONTAINS
!******************************************************************************
! Routine to execute numerical tests for debugging
!!****** fist_debug/debug_control [1.0] *
!!
!! NAME
!! debug_control
!!
!! FUNCTION
!! Routine to execute numerical tests for debugging
!!
!! AUTHOR
!! CJM
!!
!! MODIFICATION HISTORY
!! Changed Formats (JGH) 22-Nov-2000
!! Harald Forbert (Dec-2000): Changes for multiple linked lists
!!
!! SOURCE
!******************************************************************************
SUBROUTINE debug_control ( globenv, ewald_param, part, pnode, molecule, box, &
thermo, potparm, ensemble )
thermo, potparm, ensemble, lldata )
IMPLICIT NONE
@ -50,6 +84,7 @@ SUBROUTINE debug_control ( globenv, ewald_param, part, pnode, molecule, box, &
TYPE ( cell_type ), INTENT ( INOUT ) :: box
TYPE ( thermodynamic_type ), INTENT ( INOUT ) :: thermo
CHARACTER ( LEN = * ), INTENT ( IN ) :: ensemble
TYPE ( linklist_internal_data_type ) :: lldata
! Locals
TYPE ( debug_variables_type ) :: dbg
@ -92,10 +127,12 @@ SUBROUTINE debug_control ( globenv, ewald_param, part, pnode, molecule, box, &
IF ( isos /= 0 ) CALL stop_memory ( "debug_control", "f_bond", 3 * natoms )
ALLOCATE ( dbg % f_bend ( 3, natoms ), STAT = isos )
IF ( isos /= 0 ) CALL stop_memory ( "debug_control", "f_bend", 3 * natoms )
ALLOCATE ( dbg % f_torsion ( 3, natoms ), STAT = isos )
IF ( isos /= 0 ) CALL stop_memory ( "debug_control", "f_torsion", 3 * natoms )
! call to analytical nonbond force routine:
CALL force_control ( molecule, pnode, part, box, thermo, &
potparm, ewald_param, ensemble, globenv, dbg )
potparm, ewald_param, ensemble, globenv, lldata, dbg )
! Debug real-space nonbond
WRITE ( iw, '( A )' ) &
@ -135,17 +172,17 @@ SUBROUTINE debug_control ( globenv, ewald_param, part, pnode, molecule, box, &
dum = MAXLOC ( diff )
WRITE ( iw, '( A, T71, I10 )' ) &
' NON BOND PARTICLE WITH MAX ABSOLUTE ERROR IS ', dum(2)
WRITE ( iw, '( A, T21, 3G20.14 )' ) &
WRITE ( iw, '( A, T21, 3G20.12 )' ) &
' F ANAL NONBOND =', dbg % f_nonbond ( :, dum ( 2 ) )
WRITE ( iw, '( A,T21,3G20.14 )' ) &
WRITE ( iw, '( A,T21,3G20.12 )' ) &
' F NUMR NONBOND =', f_numer ( :, dum ( 2 ) )
dum = MAXLOC ( rel )
WRITE ( iw, '( A,T71,I10 )' ) &
' NON BOND PARTICLE WITH MAX RELATIVE ERROR IS ', dum(2)
WRITE ( iw, '( A,T21,3G20.14 )' ) &
WRITE ( iw, '( A,T21,3G20.12 )' ) &
' F ANAL NONBOND =', dbg % f_nonbond ( :, dum ( 2 ) )
WRITE ( iw, '( A,T21,3G20.14 )' ) &
WRITE ( iw, '( A,T21,3G20.12 )' ) &
' F NUMR NONBOND =', f_numer ( :, dum ( 2 ) )
END IF
@ -164,9 +201,9 @@ SUBROUTINE debug_control ( globenv, ewald_param, part, pnode, molecule, box, &
WRITE ( iw, '( A, T71, F10.6 )' ) ' DELTA = ', delta
END IF
CALL force_bond_numer ( molecule, f_numer, energy_numer, delta )
WRITE ( iw, '( A, T61, G20.14 )' ) ' BOND NUMER ENERGY = ', &
WRITE ( iw, '( A, T61, G20.12 )' ) ' BOND NUMER ENERGY = ', &
energy_numer
WRITE ( iw, '( A, T61, G20.14 )' ) ' BOND ANAL ENERGY = ', &
WRITE ( iw, '( A, T61, G20.12 )' ) ' BOND ANAL ENERGY = ', &
dbg % pot_bond
! computing the absolute value of the differences in the forces
@ -174,9 +211,9 @@ SUBROUTINE debug_control ( globenv, ewald_param, part, pnode, molecule, box, &
rel = diff / dbg % f_bond
! find the maximum difference and the relative and absolute errors.
WRITE ( iw, '( A,T61,G20.14 )' ) ' MAXIMUM ABSOLUTE ERROR = ', &
WRITE ( iw, '( A,T61,G20.12 )' ) ' MAXIMUM ABSOLUTE ERROR = ', &
MAXVAL ( diff )
WRITE ( iw, '( A,T61,G20.14 )' ) ' MAXIMUM RELATIVE ERROR = ', &
WRITE ( iw, '( A,T61,G20.12 )' ) ' MAXIMUM RELATIVE ERROR = ', &
MAXVAL ( rel )
! write out the particle number and forces of
@ -184,17 +221,17 @@ SUBROUTINE debug_control ( globenv, ewald_param, part, pnode, molecule, box, &
dum = MAXLOC ( diff )
WRITE ( iw, '( A,T71,I10 )' ) &
' NON BOND PARTICLE WITH MAX ABSOLUTE ERROR IS ', dum(2)
WRITE ( iw, '( A,T21,3G20.14 )' ) ' F ANAL BOND =', &
WRITE ( iw, '( A,T21,3G20.12 )' ) ' F ANAL BOND =', &
dbg % f_bond(:,dum(2))
WRITE ( iw, '( A,T21,3G20.14 )' ) ' F NUMR BOND =', &
WRITE ( iw, '( A,T21,3G20.12 )' ) ' F NUMR BOND =', &
f_numer(:,dum(2))
dum = MAXLOC ( rel )
WRITE ( iw, '( A,T71,I10 )' ) &
' NON BOND PARTICLE WITH MAX RELATIVE ERROR IS ', dum(2)
WRITE ( iw, '( A,T21,3G20.14 )' ) 'F ANAL BOND =', &
WRITE ( iw, '( A,T21,3G20.12 )' ) ' F ANAL BOND =', &
dbg % f_bond(:,dum(2))
WRITE ( iw, '( A,T21,3G20.14 )' ) 'F NUMR BOND =', &
WRITE ( iw, '( A,T21,3G20.12 )' ) ' F NUMR BOND =', &
f_numer(:,dum(2))
END IF
@ -219,11 +256,11 @@ SUBROUTINE debug_control ( globenv, ewald_param, part, pnode, molecule, box, &
! writing out the virials
WRITE ( iw, '( A )' ) ' PV NUMERICAL'
DO i = 1, 3
WRITE ( iw, '( T21, 3G20.14 )' ) pv_test ( i, : )
WRITE ( iw, '( T21, 3G20.12 )' ) pv_test ( i, : )
END DO
WRITE ( iw, '( A )' ) ' PV BOND'
DO i = 1, 3
WRITE ( iw, '( T21,3G20.14 )' ) dbg % pv_bond(i,:)
WRITE ( iw, '( T21,3G20.12 )' ) dbg % pv_bond(i,:)
END DO
END IF
END IF
@ -246,9 +283,9 @@ SUBROUTINE debug_control ( globenv, ewald_param, part, pnode, molecule, box, &
END IF
CALL force_bend_numer ( molecule, f_numer, energy_numer, delta )
WRITE ( iw, '( A, T61, G20.14 )' ) ' BEND NUMER ENERGY = ', &
WRITE ( iw, '( A, T61, G20.12 )' ) ' BEND NUMER ENERGY = ', &
energy_numer
WRITE ( iw, '( A, T61, G20.14 )' ) ' BEND ANAL ENERGY = ', &
WRITE ( iw, '( A, T61, G20.12 )' ) ' BEND ANAL ENERGY = ', &
dbg % pot_bend
! computing the absolute value of the differences in the forces
@ -256,9 +293,9 @@ SUBROUTINE debug_control ( globenv, ewald_param, part, pnode, molecule, box, &
rel = diff / dbg % f_bend
! find the maximum difference and the relative and absolute errors.
WRITE ( iw, '( A,T61,G20.14 )' ) ' MAXIMUM ABSOLUTE ERROR = ', &
WRITE ( iw, '( A,T61,G20.12 )' ) ' MAXIMUM ABSOLUTE ERROR = ', &
MAXVAL ( diff )
WRITE ( iw, '( A,T61,G20.14 )' ) ' MAXIMUM RELATIVE ERROR = ', &
WRITE ( iw, '( A,T61,G20.12 )' ) ' MAXIMUM RELATIVE ERROR = ', &
MAXVAL ( rel )
! write out the particle number and forces of
@ -266,17 +303,17 @@ SUBROUTINE debug_control ( globenv, ewald_param, part, pnode, molecule, box, &
dum = MAXLOC ( diff )
WRITE ( iw, '( A,T71,I10 )' ) &
' NON BEND PARTICLE WITH MAX ABSOLUTE ERROR IS ', dum(2)
WRITE ( iw, '( A,T21,3G20.14 )' ) 'F ANAL BEND =', &
WRITE ( iw, '( A,T21,3G20.12 )' ) ' F ANAL BEND =', &
dbg % f_bend(:,dum(2))
WRITE ( iw, '( A,T21,3G20.14 )' ) 'F NUMR BEND =', &
WRITE ( iw, '( A,T21,3G20.12 )' ) ' F NUMR BEND =', &
f_numer(:,dum(2))
dum = MAXLOC ( rel )
WRITE ( iw, '( A,T71,I10 )' ) &
' NON BEND PARTICLE WITH MAX RELATIVE ERROR IS ', dum(2)
WRITE ( iw, '( A,T21,3G20.14 )' ) 'F ANAL BEND =', &
WRITE ( iw, '( A,T21,3G20.12 )' ) ' F ANAL BEND =', &
dbg % f_bend(:,dum(2))
WRITE ( iw, '( A,T21,3G20.14 )' ) 'F NUMR BEND =', &
WRITE ( iw, '( A,T21,3G20.12 )' ) ' F NUMR BEND =', &
f_numer(:,dum(2))
END IF
@ -304,11 +341,11 @@ SUBROUTINE debug_control ( globenv, ewald_param, part, pnode, molecule, box, &
! writing out the virials
WRITE ( iw, '( A )' ) ' PV NUMERICAL'
DO i = 1, 3
WRITE ( iw, '( T21, 3G20.14 )' ) pv_test(i,:)
WRITE ( iw, '( T21, 3G20.12 )' ) pv_test(i,:)
END DO
WRITE ( iw, '( A )' ) ' PV BEND'
DO i = 1, 3
WRITE ( iw, '( T21, 3G20.14 )' ) dbg % pv_bend(i,:)
WRITE ( iw, '( T21, 3G20.12 )' ) dbg % pv_bend(i,:)
END DO
END IF
END IF
@ -332,16 +369,16 @@ SUBROUTINE debug_control ( globenv, ewald_param, part, pnode, molecule, box, &
' DO YOU WANT TO DEBUG YOUR G-SPACE ENERGIES (1=yes)?'
READ ( ir, * ) iflag
IF ( iflag == 1 ) THEN
WRITE ( iw, '( A, T71, I10 )' ) 'TOTAL NUMBER OF G-VECTORS= ', &
WRITE ( iw, '( A, T71, I10 )' ) ' TOTAL NUMBER OF G-VECTORS= ', &
ewald_grid % ngpts_cut
WRITE ( iw, '( A, T71, F10.4 )' ) 'ALPHA= ', &
WRITE ( iw, '( A, T71, F10.4 )' ) ' ALPHA= ', &
ewald_param % alpha
CALL energy_recip_numer ( ewald_param, pnode, box, &
ewald_grid, energy_numer, ewald_grid % ngpts_cut )
WRITE ( iw, '( A, T61, G20.14 )' ) 'G-SPACE ANAL ENERGY = ', &
ewald_grid, energy_numer, ewald_grid % ngpts_cut )
WRITE ( iw, '( A, T61, G20.12 )' ) ' G-SPACE ANAL ENERGY = ', &
dbg % pot_g
WRITE ( iw, '( A, T61, G20.14 )' ) &
'G-SPACE NUMERICAL ENERGY = ', energy_numer
WRITE ( iw, '( A, T61, G20.12 )' ) &
' G-SPACE NUMERICAL ENERGY = ', energy_numer
END IF
WRITE ( iw, '( A )' ) &
' DO YOU WANT TO DEBUG YOUR G-SPACE FORCES (1=yes)?'
@ -358,16 +395,16 @@ SUBROUTINE debug_control ( globenv, ewald_param, part, pnode, molecule, box, &
END IF
CALL force_recip_numer ( ewald_param, pnode, box, ewald_grid, &
delta, f_numer )
delta, f_numer )
! computing the absolute value of the differences in the forces
diff = ABS ( dbg % f_g - f_numer )
rel = diff / dbg % f_g
! find the maximum difference and the relative and absolute errors.
WRITE ( iw, '( A,T61,G20.14 )' ) &
WRITE ( iw, '( A,T61,G20.12 )' ) &
'MAXIMUM ABSOLUTE ERROR = ', maxval(diff)
WRITE ( iw, '( A,T61,G20.14 )' ) &
WRITE ( iw, '( A,T61,G20.12 )' ) &
'MAXIMUM RELATIVE ERROR = ', maxval(rel)
!
! write out the particle number and forces of the max absolute
@ -375,17 +412,17 @@ SUBROUTINE debug_control ( globenv, ewald_param, part, pnode, molecule, box, &
dum = MAXLOC ( diff )
WRITE ( iw, '( A,T71,I10 )' ) &
' PARTICLE WITH MAX ABSOLUTE ERROR IS ', dum(2)
WRITE ( iw, '( A,T21,3G20.14 )' ) ' F_ANAL G-SPACE =', &
WRITE ( iw, '( A,T21,3G20.12 )' ) ' F_ANAL G-SPACE =', &
dbg % f_g(:,dum(2))
WRITE ( iw, '( A,T21,3G20.14 )' ) ' F_NUMR G-SPACE =', &
WRITE ( iw, '( A,T21,3G20.12 )' ) ' F_NUMR G-SPACE =', &
f_numer(:,dum(2))
dum = MAXLOC ( rel )
WRITE ( iw, '( A,T71,I10 )' ) &
'PARTICLE WITH MAX RELATIVE ERROR IS ', dum(2)
WRITE ( iw, '( A,T21,3G20.14 )' ) ' F_ANAL G-SPACE =', &
WRITE ( iw, '( A,T21,3G20.12 )' ) ' F_ANAL G-SPACE =', &
dbg % f_g(:,dum(2))
WRITE ( iw, '( A,T21,3G20.14 )' ) ' F_NUMR G-SPACE =', &
WRITE ( iw, '( A,T21,3G20.12 )' ) ' F_NUMR G-SPACE =', &
f_numer(:,dum(2))
END IF
@ -405,16 +442,16 @@ SUBROUTINE debug_control ( globenv, ewald_param, part, pnode, molecule, box, &
WRITE ( iw, '( A,T71,F10.6 )' ) ' DELTA = ', delta
END IF
CALL pvg_numer ( ewald_param, pnode, box, &
ewald_grid, pv_test, delta )
ewald_grid, pv_test, delta )
! writing out the virials
WRITE ( iw, '( A )' ) ' PV G-SPACE NUMERICAL'
DO i = 1, 3
WRITE ( iw, '( T21,3G20.14 )' ) pv_test(i,:)
WRITE ( iw, '( T21,3G20.12 )' ) pv_test(i,:)
END DO
WRITE ( iw, '( A )' ) ' PV G-SPACE'
DO i = 1, 3
WRITE ( iw, '( T21,3G20.14 )' ) dbg % pv_g(i,:)
WRITE ( iw, '( T21,3G20.12 )' ) dbg % pv_g(i,:)
END DO
END IF
END IF
@ -427,7 +464,7 @@ SUBROUTINE debug_control ( globenv, ewald_param, part, pnode, molecule, box, &
IF ( iflag == 1 ) THEN
!
! get numerical virial
WRITE ( iw, '( A )' ) 'ENTER A DELTA LESS THAN 1'
WRITE ( iw, '( A )' ) ' ENTER A DELTA LESS THAN 1'
READ ( ir, * ) delta
IF ( delta >= 1.0_dbl .OR. delta <= 0.0_dbl ) THEN
delta = 1.0E-5_dbl
@ -441,11 +478,11 @@ SUBROUTINE debug_control ( globenv, ewald_param, part, pnode, molecule, box, &
! writing out the virials
WRITE ( iw, '( A )' ) ' PV NUMERICAL'
DO i = 1, 3
WRITE ( iw, '( T21,3G20.14 )' ) pv_test(i,:)
WRITE ( iw, '( T21,3G20.12 )' ) pv_test(i,:)
END DO
WRITE ( iw, '( A )' ) ' PV NONBOND'
DO i = 1, 3
WRITE ( iw, '( T21,3G20.14 )' ) dbg % pv_nonbond(i,:)
WRITE ( iw, '( T21,3G20.12 )' ) dbg % pv_nonbond(i,:)
END DO
END IF
@ -458,6 +495,7 @@ SUBROUTINE debug_control ( globenv, ewald_param, part, pnode, molecule, box, &
'***************************************************'
END SUBROUTINE debug_control
!!*****
!******************************************************************************
END MODULE fist_debug

121
src/fist_force.F
View file

@ -2,6 +2,22 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C ) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/fist_force [1.0] *
!!
!! NAME
!! fist_force
!!
!! FUNCTION
!!
!! AUTHOR
!! CJM & JGH
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE fist_force
USE dg_types, ONLY : dg_type
@ -15,7 +31,8 @@ MODULE fist_force
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_neighbor, linklist_exclusion, linklist_bonds, &
linklist_bends, linklist_torsions
USE message_passing, ONLY : mp_sum
USE pair_potential, ONLY : potentialparm_type
USE particle_types, ONLY : particle_type
@ -25,26 +42,45 @@ MODULE fist_force
USE simulation_cell, ONLY : cell_type, get_hinv
USE stop_program, ONLY : stop_prg, 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
REAL ( dbl ) :: pot_nonbond, pot_g, pot_bond, pot_bend, pot_torsion
REAL ( dbl ), DIMENSION ( :, : ), POINTER :: &
f_nonbond, f_g, f_bond, f_bend
REAL ( dbl ), DIMENSION ( 3, 3 ) :: pv_nonbond, pv_g, pv_bond, pv_bend
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
!******************************************************************************
!! Calculates the total potential energy, total force, and the
!! total pressure tensor from the potentials
!!****** fist_force/force_control [1.0] *
!!
!! NAME
!! force_control
!!
!! FUNCTION
!! Calculates the total potential energy, total force, and the
!! total pressure tensor from the potentials
!!
!! AUTHOR
!! CJM & JGH
!!
!! MODIFICATION HISTORY
!! Harald Forbert (Dec-2000): Changes for multiple linked lists
!!
!! SOURCE
!******************************************************************************
SUBROUTINE force_control ( molecule, pnode, part, box, thermo, &
potparm, ewald_param, ensemble, fc_global, debug )
potparm, ewald_param, ensemble, fc_global, lldata, debug)
IMPLICIT NONE
@ -58,14 +94,16 @@ SUBROUTINE force_control ( molecule, pnode, part, box, thermo, &
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), INTENT (INOUT) :: lldata
TYPE ( debug_variables_type ), INTENT ( OUT ), OPTIONAL :: debug
! Locals
INTEGER :: id, i, ii, natoms, nnodes, handle, isos
REAL ( dbl ) :: pot_nonbond, pot_bond, pot_bend, vg_coulomb
REAL ( dbl ) :: pot_nonbond, pot_bond, pot_bend, pot_torsion, vg_coulomb
REAL ( dbl ), DIMENSION ( :,: ), ALLOCATABLE, SAVE :: f_nonbond
REAL ( dbl ), DIMENSION ( 3,3 ) :: pv_nonbond, pv_bond, pv_bend
REAL ( dbl ), DIMENSION ( 3,3 ) :: pv_nonbond, pv_bond, pv_bend, pv_torsion
REAL ( dbl ), DIMENSION ( :,: ), ALLOCATABLE :: fg_coulomb
REAL ( dbl ), DIMENSION ( :,: ), ALLOCATABLE :: f_total
REAL ( dbl ), DIMENSION ( 3,3 ) :: pv_g, ident
REAL ( dbl ), DIMENSION ( 3,3 ) :: pv_bc
TYPE ( pw_grid_type ), SAVE :: grid_s, grid_b, grid_ewald
@ -122,7 +160,7 @@ SUBROUTINE force_control ( molecule, pnode, part, box, thermo, &
!
! first check with list_control to update neighbor lists
!
CALL list_control ( pnode, part, box )
CALL list_control ( lldata, pnode, part, box )
!
! initial force, energy and pressure tensor arrays
!
@ -133,6 +171,7 @@ SUBROUTINE force_control ( molecule, pnode, part, box, thermo, &
END DO
pv_bond = 0.0_dbl
pv_bend = 0.0_dbl
pv_torsion = 0.0_dbl
thermo % pot = 0.0_dbl
thermo % pv = 0.0_dbl
!
@ -175,23 +214,24 @@ SUBROUTINE force_control ( molecule, pnode, part, box, thermo, &
! get intramolecular forces
!
IF ( PRESENT ( debug ) ) THEN
CALL force_intra_control ( molecule, pot_bond, pot_bend, &
pv_bond, pv_bend, debug % f_bond,debug % f_bend )
CALL force_intra_control ( molecule, pot_bond, pot_bend, pot_torsion, &
pv_bond, pv_bend, pv_torsion, &
debug % f_bond, debug % f_bend, debug % f_torsion )
ELSE
CALL force_intra_control ( molecule, pot_bond, pot_bend, &
pv_bond, pv_bend )
CALL force_intra_control ( molecule, pot_bond, pot_bend, pot_torsion, &
pv_bond, pv_bend, pv_torsion )
END IF
!
! add up all the potential energies
!
IF ( ewald_param % ewald_type == 'NONE' ) THEN
thermo % pot = pot_nonbond + pot_bond + pot_bend
thermo % pot = pot_nonbond + pot_bond + pot_bend + pot_torsion
#if defined ( __parallel )
CALL mp_sum ( thermo % pot,fc_global % group )
#endif
ELSE
thermo % pot = pot_nonbond + pot_bond + pot_bend + thermo % e_bonded
thermo % pot = pot_nonbond + pot_bond + pot_bend + pot_torsion + thermo % e_bonded
thermo % gspace = vg_coulomb
#if defined ( __parallel )
CALL mp_sum ( thermo % pot,fc_global % group )
@ -208,23 +248,26 @@ SUBROUTINE force_control ( molecule, pnode, part, box, thermo, &
! 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_nonbond ( 1, i ) = part ( i ) % f ( 1 ) + f_nonbond ( 1, i )
f_nonbond ( 2, i ) = part ( i ) % f ( 2 ) + f_nonbond ( 2, i )
f_nonbond ( 3, i ) = part ( i ) % f ( 3 ) + f_nonbond ( 3, i )
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
IF ( ewald_param % ewald_type /= 'NONE' ) THEN
DO i = 1, nnodes
ii = pnode ( i ) % p % iatom
f_nonbond ( 1, ii ) = f_nonbond ( 1, ii ) + fg_coulomb ( 1, i )
f_nonbond ( 2, ii ) = f_nonbond ( 2, ii ) + fg_coulomb ( 2, i )
f_nonbond ( 3, ii ) = f_nonbond ( 3, ii ) + fg_coulomb ( 3, i )
f_total ( 1, ii ) = f_total ( 1, ii ) + fg_coulomb ( 1, i )
f_total ( 2, ii ) = f_total ( 2, ii ) + fg_coulomb ( 2, i )
f_total ( 3, ii ) = f_total ( 3, ii ) + fg_coulomb ( 3, i )
END DO
END IF
! add up all the pressure tensors
IF ( ewald_param % ewald_type == 'NONE' ) THEN
thermo % pv = pv_nonbond + pv_bond + pv_bend
thermo % pv = pv_nonbond + pv_bond + pv_bend + pv_torsion
#if defined ( __parallel )
CALL mp_sum ( thermo % pv, fc_global % group )
#endif
@ -234,7 +277,7 @@ SUBROUTINE force_control ( molecule, pnode, part, box, thermo, &
ident ( i, i ) = 1.0_dbl
END DO
thermo % pv = pv_nonbond + pv_bond + pv_bend + pv_bc
thermo % pv = pv_nonbond + pv_bond + pv_bend + pv_torsion + pv_bc
#if defined ( __parallel )
CALL mp_sum ( thermo % pv,fc_global % group )
#endif
@ -249,11 +292,13 @@ SUBROUTINE force_control ( molecule, pnode, part, box, thermo, &
IF ( PRESENT ( debug ) ) THEN
debug % pot_bond = pot_bond
debug % pot_bend = pot_bend
debug % pot_torsion = pot_torsion
debug % pot_nonbond = pot_nonbond
debug % f_nonbond = f_nonbond
debug % pv_nonbond = pv_nonbond
debug % pv_bond = pv_bond
debug % pv_bend = pv_bend
debug % pv_torsion = pv_torsion
IF ( ewald_param % ewald_type /= 'NONE' ) THEN
debug % pot_g = vg_coulomb
debug % f_g = fg_coulomb
@ -265,25 +310,31 @@ SUBROUTINE force_control ( molecule, pnode, part, box, thermo, &
END IF
END IF
#if defined ( __parallel )
CALL mp_sum ( f_total, fc_global % group )
#endif
DO i = 1, natoms
part ( i ) % f ( 1 ) = f_total ( 1, i )
part ( i ) % f ( 2 ) = f_total ( 2, i )
part ( i ) % f ( 3 ) = f_total ( 3, i )
END DO
! deallocating all local variables
isos = 0
IF ( ALLOCATED ( fg_coulomb ) ) DEALLOCATE ( fg_coulomb, STAT = isos )
IF ( isos /= 0 ) CALL stop_memory ( 'force_control', 'fg_coulomb' )
#if defined ( __parallel )
CALL mp_sum ( f_nonbond, fc_global % group )
#endif
isos = 0
IF ( ALLOCATED ( f_total ) ) DEALLOCATE ( f_total, STAT = isos )
IF ( isos /= 0 ) CALL stop_memory ( 'force_control', 'f_total' )
DO i = 1, natoms
part ( i ) % f ( 1 ) = f_nonbond ( 1, i )
part ( i ) % f ( 2 ) = f_nonbond ( 2, i )
part ( i ) % f ( 3 ) = f_nonbond ( 3, i )
END DO
CALL timestop ( zero, handle )
END SUBROUTINE force_control
!!*****
!******************************************************************************
END MODULE fist_force
!******************************************************************************

70
src/fist_force_numer.F
View file

@ -2,6 +2,21 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/fist_force_numer [1.0] *
!!
!! NAME
!! fist_force_numer
!!
!! FUNCTION
!!
!! AUTHOR
!! CJM
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE fist_force_numer
@ -23,6 +38,9 @@ MODULE fist_force_numer
PUBLIC :: force_recip_numer, pvbond_numer, pvbend_numer, ptens_numer
PUBLIC :: pvg_numer, de_g_numer, potential_g_numer, energy_recip_numer
!!*****
!******************************************************************************
CONTAINS
!******************************************************************************
@ -168,7 +186,7 @@ SUBROUTINE force_bend_numer ( molecule, f_bend_numer, v_bend_numer, &
DO idim = 1, 3
delta(idim) = numerical_shift
b32_minus_delta = b32 - delta
d32m = sqrt(dot_product(b32_minus_delta,b12_minus_delta))
d32m = sqrt(dot_product(b32_minus_delta,b32_minus_delta))
id32m = 1.0_dbl/d32m
dist = dot_product(b32_minus_delta,b12)
theta = acos(dist*id32m*id12)
@ -300,7 +318,7 @@ SUBROUTINE force_nonbond_numer ( ewald_param, pnode, box, potparm, &
delta(id) = numerical_shift
rij_minus_delta_sq = dot_product(rij-delta,rij-delta)
rij_plus_delta_sq = dot_product(rij+delta,rij+delta)
IF (qi==0.AND.qj==0) THEN
IF (qi==0.AND.qj==0) THEN
CALL potential_f(rij_minus_delta_sq,potparm,qi,qj, &
iatomtype, jatomtype, energy_minus )
CALL potential_f(rij_plus_delta_sq,potparm,qi,qj, &
@ -315,13 +333,13 @@ SUBROUTINE force_nonbond_numer ( ewald_param, pnode, box, potparm, &
f_numer(id,j) = f_numer(id,j) - (energy_plus-energy_minus)
delta = 0.0_dbl
END DO DIM_LOOP
IF (qi==0.AND.qj==0) THEN
IF (qi==0.AND.qj==0) THEN
CALL potential_f ( rijsq, potparm, qi, qj, &
iatomtype, jatomtype, energy )
ELSE
CALL potential_f ( rijsq, potparm, qi, qj, &
iatomtype, jatomtype, energy, ewald_param )
ENDIF
ENDIF
energy_numer = energy_numer + energy
END IF NOTMATCH
@ -353,17 +371,17 @@ SUBROUTINE force_nonbond_numer ( ewald_param, pnode, box, potparm, &
delta(id) = numerical_shift
rij_minus_delta_sq = dot_product(rij-delta,rij-delta)
rij_plus_delta_sq = dot_product(rij+delta,rij+delta)
IF (qi==0.AND.qj==0) THEN
IF (qi==0.AND.qj==0) THEN
CALL potential_f(rij_minus_delta_sq,potparm,qi,qj, &
iatomtype, jatomtype, energy_minus )
CALL potential_f(rij_plus_delta_sq,potparm,qi,qj, &
iatomtype, jatomtype, energy_plus )
ELSE
ELSE
CALL potential_f(rij_minus_delta_sq,potparm,qi,qj, &
iatomtype, jatomtype, energy_minus, ewald_param )
CALL potential_f(rij_plus_delta_sq,potparm,qi,qj, &
iatomtype, jatomtype, energy_plus, ewald_param )
ENDIF
ENDIF
f_numer(id,i) = f_numer(id,i) &
+ (energy_plus-energy_minus)
f_numer(id,j) = f_numer(id,j) &
@ -371,13 +389,13 @@ SUBROUTINE force_nonbond_numer ( ewald_param, pnode, box, potparm, &
delta = 0.0_dbl
END DO DIM_LOOP2
IF (qi==0.AND.qj==0) THEN
IF (qi==0.AND.qj==0) THEN
CALL potential_f ( rijsq, potparm, qi, qj, &
iatomtype, jatomtype, energy )
ELSE
CALL potential_f ( rijsq, potparm, qi, qj, &
iatomtype, jatomtype, energy, ewald_param )
ENDIF
ENDIF
energy_numer = energy_numer + energy
END IF
END DO
@ -403,17 +421,17 @@ SUBROUTINE force_nonbond_numer ( ewald_param, pnode, box, potparm, &
delta(id) = numerical_shift
rij_minus_delta_sq = dot_product(rij-delta,rij-delta)
rij_plus_delta_sq = dot_product(rij+delta,rij+delta)
IF (qi==0.AND.qj==0) THEN
IF (qi==0.AND.qj==0) THEN
CALL potential_f(rij_minus_delta_sq,potparm,qi,qj, &
iatomtype, jatomtype, energy_minus )
CALL potential_f(rij_plus_delta_sq,potparm,qi,qj, &
iatomtype, jatomtype, energy_plus )
ELSE
ELSE
CALL potential_f(rij_minus_delta_sq,potparm,qi,qj, &
iatomtype, jatomtype, energy_minus, ewald_param )
CALL potential_f(rij_plus_delta_sq,potparm,qi,qj, &
iatomtype, jatomtype, energy_plus, ewald_param )
ENDIF
ENDIF
f_numer(id,i) = f_numer(id,i) &
+ (energy_plus-energy_minus)
f_numer(id,i) = f_numer(id,i) &
@ -421,13 +439,13 @@ SUBROUTINE force_nonbond_numer ( ewald_param, pnode, box, potparm, &
delta = 0.0_dbl
END DO DIM_LOOP3
IF (qi==0.AND.qj==0) THEN
IF (qi==0.AND.qj==0) THEN
CALL potential_f ( rijsq, potparm, qi, qj, &
iatomtype, jatomtype, energy )
ELSE
CALL potential_f ( rijsq, potparm, qi, qj, &
iatomtype, jatomtype, energy, ewald_param )
ENDIF
ENDIF
energy_numer = energy_numer + energy
END IF
END DO
@ -498,10 +516,10 @@ SUBROUTINE force_recip_numer ( ewald_param, pnode, box, pw_grid, &
DO idim = 1, 3
r_delta(idim) = r_delta(idim) + numerical_shift
CALL de_g_numer(ep, sum_igr, gauss, pnode(i) % p % r , &
r_delta, charge, pw_grid % g, ngtot)
r_delta, charge, pw_grid % g, ngtot)
r_delta(idim) = r_delta(idim) - 2.0_dbl*numerical_shift
CALL de_g_numer(em, sum_igr, gauss, pnode(i) % p % r , &
r_delta, charge, pw_grid % g, ngtot)
r_delta, charge, pw_grid % g, ngtot)
f_numer(idim,i) = ( em - ep ) / epsilon0 / box % deth
r_delta(idim) = r_delta(idim) + numerical_shift
END DO
@ -1217,13 +1235,13 @@ SUBROUTINE getv_ptens ( ewald_param, pnode, x, vnb, box, potparm )
END DO EXCL
NOTMATCH: IF ( .NOT. match) THEN
IF (qi==0.AND.qj==0) THEN
IF (qi==0.AND.qj==0) THEN
CALL potential_f ( rijsq, potparm, qi, qj, &
iatomtype, jatomtype, energy )
ELSE
ELSE
CALL potential_f ( rijsq, potparm, qi, qj, &
iatomtype, jatomtype, energy, ewald_param )
ENDIF
ENDIF
!
! summing up the potential energy
!
@ -1252,13 +1270,13 @@ SUBROUTINE getv_ptens ( ewald_param, pnode, x, vnb, box, potparm )
CALL find_image(s,perd,vec,box % hmat,rijsq,rij)
IF ( rijsq <= potparm ( iatomtype, jatomtype ) % rcutsq ) &
THEN
IF (qi==0.AND.qj==0) THEN
IF (qi==0.AND.qj==0) THEN
CALL potential_f ( rijsq, potparm, qi, qj, &
iatomtype, jatomtype, energy )
ELSE
ELSE
CALL potential_f ( rijsq, potparm, qi, qj, &
iatomtype, jatomtype, energy, ewald_param )
ENDIF
ENDIF
vnb = vnb + energy
END IF
END DO
@ -1280,13 +1298,13 @@ SUBROUTINE getv_ptens ( ewald_param, pnode, x, vnb, box, potparm )
CALL find_image ( s, perd, vec, box % hmat, rijsq, rij )
IF ( rijsq <= potparm ( iatomtype, iatomtype ) % rcutsq ) &
THEN
IF (qi==0.AND.qj==0) THEN
IF (qi==0.AND.qj==0) THEN
CALL potential_f ( rijsq, potparm, qi, qj, &
iatomtype, jatomtype, energy )
ELSE
ELSE
CALL potential_f ( rijsq, potparm, qi, qj, &
iatomtype, jatomtype, energy, ewald_param )
ENDIF
ENDIF
vnb = vnb + energy
END IF
@ -1303,3 +1321,5 @@ END SUBROUTINE getv_ptens
!******************************************************************************
END MODULE fist_force_numer
!******************************************************************************

39
src/fist_global.F
View file

@ -1,39 +0,0 @@
!-----------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!! Definition of some global variables for the FIST code
MODULE fist_global
USE global_types, ONLY : global_environment_type
USE kinds, ONLY : dbl
IMPLICIT NONE
PRIVATE
PUBLIC :: set_fist_global, fistpar
TYPE ( global_environment_type ) :: fistpar
CONTAINS
!******************************************************************************
SUBROUTINE set_fist_global ( glopar )
IMPLICIT NONE
! Arguments
TYPE ( global_environment_type ), INTENT ( IN ) :: glopar
!------------------------------------------------------------------------------
fistpar = glopar
END SUBROUTINE set_fist_global
!******************************************************************************
END MODULE fist_global

65
src/fist_input.F
View file

@ -2,6 +2,21 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/fist_input [1.0] *
!!
!! NAME
!! fist_input
!!
!! FUNCTION
!!
!! AUTHOR
!! CJM & JGH
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE fist_input
@ -16,22 +31,43 @@ MODULE fist_input
PRIVATE
PUBLIC :: read_fist_section
!!*****
!******************************************************************************
CONTAINS
!!>---------------------------------------------------------------------------!
!! SECTION: &fist ... &end !
!! !
!! simulation [md,debug] !
!! printlevel fistpar%print_level !
!! units [kelvin,atomic] !
!! periodic [0,1][0,1][0,1] !
!! Ewald_type [pme_gauss,ewald_gauss] !
!! Ewald_param alpha,[gmax,ns_max,epsilon] !
!! set_file "filename" !
!! input_file "filename" !
!! !
!!<---------------------------------------------------------------------------!
!******************************************************************************
!!****** fist_input/read_fist_section [1.0] *
!!
!! NAME
!! read_fist_section
!!
!! FUNCTION
!!
!! AUTHOR
!! CJM & JGH
!!
!! MODIFICATION HISTORY
!! none
!!
!! INPUTS
!!I---------------------------------------------------------------------------I
!!I SECTION: &fist ... &end I
!!I I
!!I simulation [md,debug] I
!!I printlevel fistpar%print_level I
!!I units [kelvin,atomic] I
!!I periodic [0,1][0,1][0,1] I
!!I Ewald_type [pme_gauss,ewald_gauss] I
!!I Ewald_param alpha,[gmax,ns_max,epsilon] I
!!I set_file "filename" I
!!I input_file "filename" I
!!I I
!!I---------------------------------------------------------------------------I
!!
!! SOURCE
!******************************************************************************
SUBROUTINE read_fist_section ( setup, ewald_param, fistpar )
@ -162,6 +198,9 @@ SUBROUTINE read_fist_section ( setup, ewald_param, fistpar )
END SUBROUTINE read_fist_section
!!*****
!******************************************************************************
END MODULE fist_input
!******************************************************************************

117
src/fist_intra_force.F
View file

@ -2,6 +2,21 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/fist_intra_force [1.0] *
!!
!! NAME
!! fist_intra_force
!!
!! FUNCTION
!!
!! AUTHOR
!! CJM
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE fist_intra_force
@ -9,20 +24,38 @@ MODULE fist_intra_force
USE mol_force, ONLY : force_bonds, force_bends, force_torsions, &
get_pv_bond, get_pv_bend, get_pv_torsion
USE molecule_types, ONLY : molecule_structure_type, linklist_bonds, &
linklist_bends
linklist_bends, linklist_torsions
USE timings, ONLY : timeset, timestop
IMPLICIT NONE
PRIVATE
PUBLIC :: force_intra_control
!!*****
!******************************************************************************
CONTAINS
!******************************************************************************
!!****** fist_intra_force/force_intra_control [1.0] *
!!
!! NAME
!! force_intra_control
!!
!! FUNCTION
!!
!! AUTHOR
!! CJM
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
SUBROUTINE force_intra_control ( molecule, v_bond, v_bend, pvbond, pvbend, &
f_bond, f_bend )
SUBROUTINE force_intra_control ( molecule, pot_bond, pot_bend, pot_torsion, &
pv_bond, pv_bend, pv_torsion, f_bond, f_bend, f_torsion )
! parses up intramolecular properties
@ -30,17 +63,22 @@ SUBROUTINE force_intra_control ( molecule, v_bond, v_bend, pvbond, pvbend, &
! Arguments
TYPE (molecule_structure_type ), DIMENSION ( : ), INTENT ( IN ) :: molecule
REAL ( dbl ), INTENT ( INOUT ) :: v_bond, v_bend
REAL ( dbl ), DIMENSION ( :, : ), INTENT ( INOUT ) :: pvbond, pvbend
REAL ( dbl ), DIMENSION ( :, : ), OPTIONAL, INTENT ( OUT ) :: f_bond, f_bend
REAL ( dbl ), INTENT ( INOUT ) :: pot_bond, pot_bend, pot_torsion
REAL ( dbl ), DIMENSION ( :, : ), INTENT ( INOUT ) :: &
pv_bond, pv_bend, pv_torsion
REAL ( dbl ), DIMENSION ( :, : ), OPTIONAL, INTENT ( OUT ) :: &
f_bond, f_bend, f_torsion
! Locals
REAL ( dbl ) :: d12, d32, id12, id32, dist, theta
REAL ( dbl ) :: d12, d32, id12, id32, s12, s32, is12, is32, &
dist, dist1, dist2, w1, w2, iw1, iw2, theta
REAL ( dbl ) :: energy, fscalar
REAL ( dbl ), DIMENSION (3) :: rij, b12, b32, g1, g2, g3
REAL ( dbl ), DIMENSION (3) :: rij, b12, b32, t12, t32, t43, &
u1, u2, g1, g2, g3, gt1, gt2, gt3, gt4
TYPE (linklist_bonds), POINTER :: llbond
TYPE (linklist_bends), POINTER :: llbend
INTEGER :: ibond, ibend, imol, handle
TYPE (linklist_torsions), POINTER :: lltorsion
INTEGER :: ibond, ibend, itorsion, imol, handle
!------------------------------------------------------------------------------
@ -48,24 +86,27 @@ SUBROUTINE force_intra_control ( molecule, v_bond, v_bend, pvbond, pvbend, &
IF ( PRESENT ( f_bond)) f_bond = 0._dbl
IF ( PRESENT ( f_bend)) f_bend = 0._dbl
v_bond = 0._dbl
v_bend = 0._dbl
IF ( PRESENT ( f_torsion)) f_torsion = 0._dbl
pot_bond = 0._dbl
pot_bend = 0._dbl
pot_torsion = 0._dbl
MOL: DO imol = 1, size(molecule)
llbond => molecule(imol) %ll_bonds
llbend => molecule(imol) %ll_bends
lltorsion => molecule(imol) %ll_torsions
BOND: DO ibond = 1, molecule(imol) %nbonds_mol
rij = llbond%p1%r - llbond%p2%r
CALL force_bonds(rij,llbond%bond_param%r0,llbond%bond_param%k, &
energy,fscalar)
v_bond = v_bond + energy
pot_bond = pot_bond + energy
llbond%p1%f = llbond%p1%f - rij*fscalar
llbond%p2%f = llbond%p2%f + rij*fscalar
! computing the pressure tensor
CALL get_pv_bond ( -rij * fscalar, rij * fscalar, &
llbond%p1%r, llbond%p2%r, pvbond )
llbond%p1%r, llbond%p2%r, pv_bond )
! the contribution from the bonds. ONLY FOR DEBUG
IF ( PRESENT ( f_bond)) THEN
@ -90,14 +131,14 @@ SUBROUTINE force_intra_control ( molecule, v_bond, v_bend, pvbond, pvbend, &
CALL force_bends(b12,b32,d12,d32,id12,id32,dist,theta, &
llbend%bend_param%theta0,llbend%bend_param%k,g1,g2,g3,energy, &
fscalar)
v_bend = v_bend + energy
pot_bend = pot_bend + energy
llbend%p1%f = llbend%p1%f + fscalar*g1
llbend%p2%f = llbend%p2%f + fscalar*g2
llbend%p3%f = llbend%p3%f + fscalar*g3
! computing the pressure tensor
CALL get_pv_bend(fscalar*g1,fscalar*g2,fscalar*g3,llbend%p1%r, &
llbend%p2%r,llbend%p3%r,pvbend)
llbend%p2%r,llbend%p3%r,pv_bend)
! the contribution from the bends. ONLY FOR DEBUG
IF ( PRESENT ( f_bend)) THEN
@ -112,12 +153,58 @@ SUBROUTINE force_intra_control ( molecule, v_bond, v_bend, pvbond, pvbend, &
llbend => llbend%next
END DO BEND
TORSION: DO itorsion = 1, molecule(imol) %ntorsions_mol
t12 = lltorsion%p1%r - lltorsion%p2%r
t32 = lltorsion%p3%r - lltorsion%p2%r
t43 = lltorsion%p4%r - lltorsion%p3%r
s12 = sqrt(dot_product(t12,t12))
is12 = 1._dbl/s12
s32 = sqrt(dot_product(t32,t32))
is32 = 1._dbl/s32
dist1 = dot_product(t12,t32)
dist2 = dot_product(t43,t32)
u1 = t12 - dist1*t32*is32**2
u2 = t43 - dist2*t32*is32**2
w1 = sqrt(dot_product(u1,u1))
w2 = sqrt(dot_product(u2,u2))
iw1 = 1._dbl/w1
iw2 = 1._dbl/w2
CALL force_torsions(is32,dist1,dist2,u1,u2,iw1,iw2, &
lltorsion%torsion_param%a,gt1,gt2,gt3,gt4,energy,fscalar)
pot_torsion = pot_torsion + energy
lltorsion%p1%f = lltorsion%p1%f + fscalar*gt1
lltorsion%p2%f = lltorsion%p2%f + fscalar*gt2
lltorsion%p3%f = lltorsion%p3%f + fscalar*gt3
lltorsion%p4%f = lltorsion%p4%f + fscalar*gt4
! computing the pressure tensor
CALL get_pv_torsion(fscalar*gt1,fscalar*gt2,fscalar*gt3,fscalar*gt4, &
lltorsion%p1%r,lltorsion%p2%r,lltorsion%p3%r,lltorsion%p4%r,pv_torsion)
! the contribution from the torsions. ONLY FOR DEBUG
IF ( PRESENT ( f_torsion)) THEN
f_torsion(:,lltorsion%index(1)) = f_torsion(:,lltorsion%index(1)) + &
fscalar*gt1
f_torsion(:,lltorsion%index(2)) = f_torsion(:,lltorsion%index(2)) + &
fscalar*gt2
f_torsion(:,lltorsion%index(3)) = f_torsion(:,lltorsion%index(3)) + &
fscalar*gt3
f_torsion(:,lltorsion%index(4)) = f_torsion(:,lltorsion%index(4)) + &
fscalar*gt3
END IF
lltorsion => lltorsion%next
END DO TORSION
END DO MOL
CALL timestop ( 0._dbl, handle )
END SUBROUTINE force_intra_control
!!*****
!******************************************************************************
END MODULE fist_intra_force
!******************************************************************************

18
src/fist_nonbond_force.F
View file

@ -2,6 +2,21 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/fist_nonbond_force [1.0] *
!!
!! NAME
!! fist_nonbond_force
!!
!! FUNCTION
!!
!! AUTHOR
!! CJM
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE fist_nonbond_force
@ -20,6 +35,9 @@ MODULE fist_nonbond_force
PRIVATE
PUBLIC :: force_nonbond, bonded_correct_gaussian, find_image
!!*****
!******************************************************************************
CONTAINS

42
src/force_control.F
View file

@ -2,6 +2,22 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/force_control [1.0] *
!!
!! NAME
!! force_control
!!
!! FUNCTION
!! Interface routine for the force calculations
!!
!! AUTHOR
!! CJM & JGH
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE force_control
@ -18,9 +34,28 @@ MODULE force_control
PRIVATE
PUBLIC :: force
!!*****
!******************************************************************************
CONTAINS
!******************************************************************************
!!****** force_control/force [1.0] *
!!
!! NAME
!! force
!!
!! FUNCTION
!! Interface routine for force calculations
!!
!! AUTHOR
!! CJM & JGH
!!
!! MODIFICATION HISTORY
!! Harald Forbert (Dec-2000): Changes for multiple linked lists
!!
!! SOURCE
!******************************************************************************
SUBROUTINE force ( struc, inter, thermo, simpar, ewald_param, globenv )
@ -44,17 +79,20 @@ SUBROUTINE force ( struc, inter, thermo, simpar, ewald_param, globenv )
CASE ( 'FIST' )
CALL fist_force_control ( struc % molecule, struc % pnode, struc % part, &
struc % box, thermo, inter%potparm, ewald_param, simpar % ensemble, &
globenv )
globenv, struc%ll_data )
CASE ( 'TBMD' )
CALL tbmd_force_control ( struc % molecule, struc % pnode, struc % part, &
struc % box, thermo, inter%potparm, ewald_param, simpar % ensemble, &
globenv )
globenv , struc%ll_data )
END SELECT
END SUBROUTINE force
!!*****
!******************************************************************************
END MODULE force_control
!******************************************************************************

108
src/force_fields.F
View file

@ -2,6 +2,21 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/force_fields [1.0] *
!!
!! NAME
!! force_fields
!!
!! FUNCTION
!!
!! AUTHOR
!! CJM
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE force_fields
@ -32,6 +47,9 @@ MODULE force_fields
INTEGER :: scr
TYPE ( global_environment_type ) :: globenv
!!*****
!******************************************************************************
CONTAINS
@ -297,8 +315,8 @@ SUBROUTINE input_nonbonded ( potparm, atom_names, set_fn )
DO i1 = 1, natom_types
DO i2 = i1, natom_types
IF (potparm(i1,i2)%type=='LENNARD-JONES') THEN
WRITE ( iw, '( A,A,T35,A,T45,A,T58,A,T71,F10.4 )' ) &
' FORCE FIELD| ', 'Lennard-Jones(12-6)', &
WRITE ( iw, '( A,A,T30,A,T37,A,T58,A,T71,F10.4 )' ) &
' FORCE FIELD| ', 'LJ (12-6)', &
adjustl(atom_names(i1)), adjustl(atom_names(i2)), &
' epsilon [K]', potparm(i1,i2) %lj%epsilon
WRITE ( iw, '( A,T58,A,T71,F10.4 )' ) ' FORCE FIELD| ', &
@ -306,7 +324,7 @@ SUBROUTINE input_nonbonded ( potparm, atom_names, set_fn )
WRITE ( iw, '( A,T58,A,T71,F10.4 )' ) ' FORCE FIELD| ', &
' rcut [A]', sqrt(potparm(i1,i2)%rcutsq)
ELSE IF (potparm(i1,i2)%type=='WILLIAMS') THEN
WRITE ( iw, '( A,A,T30,A,T40,A,T55,A,T68,E13.4 )' ) &
WRITE ( iw, '( A,A,T30,A,T37,A,T55,A,T68,E13.4 )' ) &
' FORCE FIELD| ', 'Williams', adjustl(atom_names(i1)), &
adjustl(atom_names(i2)), 'A[K]', potparm(i1,i2) %willis%a
WRITE ( iw, '( A,T47,A,T71,F10.4 )' ) ' FORCE FIELD| ', &
@ -316,7 +334,7 @@ SUBROUTINE input_nonbonded ( potparm, atom_names, set_fn )
WRITE ( iw, '( A,T47,A,T71,F10.4 )' ) ' FORCE FIELD| ', &
' rcut [A]', sqrt(potparm(i1,i2)%rcutsq)
ELSE IF (potparm(i1,i2)%type=='GOODWIN') THEN
WRITE ( iw, '( A,A,T30,A,T40,A,T55,A,T68,E13.4 )' ) &
WRITE ( iw, '( A,A,T30,A,T37,A,T55,A,T68,E13.4 )' ) &
' FORCE FIELD| ', 'Goodwin', adjustl(atom_names(i1)), &
adjustl(atom_names(i2)), 'Ve0[eV]', &
potparm(i1,i2) %goodwin%vr0
@ -404,6 +422,7 @@ SUBROUTINE input_bonds(bond_parm,atom_names,set_fn)
CALL read_line
ilen = 0
CALL cfield(string2,ilen)
CALL uppercase(string2)
SELECT CASE (string2)
CASE DEFAULT
CALL p_error()
@ -443,12 +462,12 @@ SUBROUTINE input_bonds(bond_parm,atom_names,set_fn)
DO i1 = 1, natom_types
DO i2 = i1, natom_types
IF (index(bond_parm(i1,i2)%type,'HARMONIC') /= 0 ) THEN
WRITE ( iw, '( A, A, T30, A, T40, A, T57, F10.4, F14.2 )' ) &
' FIST| ', ADJUSTL ( bond_parm ( i1, i2 ) % type ), &
ADJUSTL ( atom_names ( i1 ) ), &
ADJUSTL ( atom_names ( i2 ) ), &
bond_parm ( i1, i2 ) % r0, &
bond_parm ( i1, i2 ) % k
WRITE ( iw, '( A, T30, A, T37, A, T58, A, T71, F10.2 )' ) &
' FORCE FIELD| harmonic', &
ADJUSTL ( atom_names ( i1 ) ), ADJUSTL ( atom_names ( i2 ) ), &
' r0 [A]', bond_parm ( i1, i2 ) % r0
WRITE ( iw, '( A, T58, A, T71, F10.2 )' ) ' FORCE FIELD| ', &
' k [K*A^-2)]', bond_parm ( i1, i2 ) % k
END IF
END DO
END DO
@ -495,11 +514,11 @@ SUBROUTINE input_bends ( bend_parm, atom_names, set_fn )
natom_types=SIZE(atom_names)
nb = 0
torad = pi / 180.0_dbl
todeg = 180.0_dbl/pi
ALLOCATE (bend_parm(natom_types,natom_types,natom_types),STAT=ios)
IF ( ios /= 0 ) CALL stop_memory ( 'input_bends', &
'bend_parm', natom_types ** 3 )
torad = pi / 180.0_dbl
todeg = 180.0_dbl/pi
bend_parm(:,:,:) %type = 'null'
bend_parm(:,:,:) %theta0 = 0._dbl
bend_parm(:,:,:) %k = 0._dbl
@ -525,6 +544,7 @@ SUBROUTINE input_bends ( bend_parm, atom_names, set_fn )
CALL read_line
ilen = 0
CALL cfield(string2,ilen)
CALL uppercase(string2)
SELECT CASE (string2)
CASE DEFAULT
CALL p_error()
@ -544,7 +564,7 @@ SUBROUTINE input_bends ( bend_parm, atom_names, set_fn )
i2 = str_search(atom_names,natom_types,at2)
IF (i2==0) CALL stop_parser( 'input_bends','atom not found')
i3 = str_search(atom_names,natom_types,at3)
IF (i2==0) CALL stop_parser( 'input_bends','atom not found')
IF (i3==0) CALL stop_parser( 'input_bends','atom not found')
bend_parm(i1,i2,i3) %type = string2
bend_parm(i1,i2,i3) %theta0 = bp1*torad
bend_parm(i1,i2,i3) %k = bp2
@ -569,11 +589,12 @@ SUBROUTINE input_bends ( bend_parm, atom_names, set_fn )
DO i2 = 1, natom_types
DO i3 = i1, natom_types
IF (index(bend_parm(i1,i2,i3)%type,'HARMONIC') /= 0 ) THEN
WRITE ( iw, '( A,A,T20,A,T30,A,T40,A,T57,F10.4,F14.2 )' ) &
' FIST| ', adjustl(bend_parm(i1,i2,i3)%type ), &
adjustl(atom_names(i1)), adjustl(atom_names(i2)), &
adjustl(atom_names(i3)), bend_parm(i1,i2,i3) %theta0* &
todeg, bend_parm(i1,i2,i3) %k
WRITE ( iw, '( A,T30,A,T37,A,T44,A,T58,A,T71,F10.2 )' ) &
' FORCE FIELD| harmonic', &
adjustl(atom_names(i1)), adjustl(atom_names(i2)), adjustl(atom_names(i3)), &
'theta0 [deg]', bend_parm(i1,i2,i3) %theta0*todeg
WRITE ( iw, '( A,T58,A,T71,F10.2 )' ) ' FORCE FIELD| ', &
' k [K]', bend_parm(i1,i2,i3) %k
END IF
END DO
END DO
@ -596,9 +617,9 @@ END SUBROUTINE input_bends
!!>---------------------------------------------------------------------------!
!! SECTION: &force_field ... &end !
!! !
!! torsion !
!! power atom1 atom2 atom3 atom4 parm[] !
!! end torsion !
!! torsions !
!! jorgensen atom1 atom2 atom3 atom4 tp[1-4] !
!! end torsions !
!! !
!!<---------------------------------------------------------------------------!
SUBROUTINE input_torsions ( torsion_parm, atom_names, set_fn )
@ -611,17 +632,15 @@ SUBROUTINE input_torsions ( torsion_parm, atom_names, set_fn )
CHARACTER ( LEN = * ), INTENT ( IN ) :: set_fn
! Locals
REAL ( dbl ) :: parm(0:6), torad, todeg
INTEGER :: i, i1, i2, i3, i4, ierror, ilen, ios
INTEGER :: nb, iw, allgrp, source, natom_types
CHARACTER ( LEN = 20 ) :: type, at1, at2, at3, at4
CHARACTER ( LEN = 40 ) :: string, string2
CHARACTER ( LEN = 12 ) :: label
REAL ( dbl ), DIMENSION (0:3) :: tp
!------------------------------------------------------------------------------
torad = pi/180._dbl
todeg = 180._dbl/pi
natom_types=SIZE(atom_names)
nb = 0
ALLOCATE ( &
@ -630,8 +649,8 @@ SUBROUTINE input_torsions ( torsion_parm, atom_names, set_fn )
IF ( ios /= 0 ) CALL stop_memory ( 'input_torsions', &
'torsion_parm', natom_types ** 4 )
torsion_parm ( :, :, :, : ) % type = 'null'
DO i = 0, SIZE ( torsion_parm(1,1,1,1) % parm ) - 1
torsion_parm(:,:,:,:) %parm(i) = 0._dbl
DO i = 0, SIZE ( torsion_parm(1,1,1,1) % a ) - 1
torsion_parm(:,:,:,:) %a(i) = 0._dbl
END DO
iw = scr
@ -650,19 +669,20 @@ SUBROUTINE input_torsions ( torsion_parm, atom_names, set_fn )
ELSE
CALL read_line
DO WHILE (test_next()/='X')
ilen = 7
ilen = 8
CALL cfield(string,ilen)
CALL uppercase ( string )
IF (index(string,'TORSION') /= 0 ) THEN
IF (index(string,'TORSIONS') /= 0 ) THEN
DO
CALL read_line
ilen = 0
CALL cfield(string2,ilen)
CALL uppercase(string2)
SELECT CASE (string2)
CASE DEFAULT
CALL p_error()
CALL stop_parser( 'input_torsion','unknown option')
CASE ( 'POWER')
CASE ( 'JORGENSEN')
nb = nb + 1
ilen = 0
CALL cfield(at1,ilen)
@ -672,8 +692,8 @@ SUBROUTINE input_torsions ( torsion_parm, atom_names, set_fn )
CALL cfield(at3,ilen)
ilen = 0
CALL cfield(at4,ilen)
DO i = 1, size(parm)
parm(i) = get_real()
DO i = 0, SIZE(tp) - 1
tp(i) = get_real()
END DO
i1 = str_search(atom_names,natom_types,at1)
IF (i1==0) CALL stop_parser( 'input_torsions', &
@ -687,8 +707,10 @@ SUBROUTINE input_torsions ( torsion_parm, atom_names, set_fn )
i4 = str_search(atom_names,natom_types,at4)
IF (i4==0) CALL stop_parser( 'input_torsions', &
'atom not found')
torsion_parm(i1,i2,i3,i4) %type = type
torsion_parm(i1,i2,i3,i4) %parm = parm
torsion_parm(i1,i2,i3,i4) %type = string2
DO i = 1, SIZE(tp) - 1
torsion_parm(i1,i2,i3,i4) %a(i) = tp(i)
END DO
torsion_parm(i4,i3,i2,i1) = torsion_parm(i1,i2,i3,i4)
CASE ( 'END')
ilen = 0
@ -706,22 +728,26 @@ SUBROUTINE input_torsions ( torsion_parm, atom_names, set_fn )
IF ( globenv % ionode ) THEN
IF ( globenv % print_level >= 0 .AND. nb > 0 ) THEN
WRITE ( iw, '( A )' ) &
' FORCE FIELD| Torsion inatom_typeseraction parameters'
' FORCE FIELD| Torsion interaction parameters'
DO i1 = 1, natom_types
DO i2 = 1, natom_types
DO i3 = i2, natom_types
DO i4 = i1, natom_types
IF ( index ( torsion_parm(i1,i2,i3,i4) % type, &
'POWER' ) == 0 ) THEN
WRITE ( iw, &
'( 2A, T20, A, T30, A, T40, A, T50, A, 7F10.4 )' )&
' FIST| ', &
ADJUSTL ( torsion_parm ( i1,i2,i3,i4 ) % type ), &
'JORGENSEN' ) /= 0 ) THEN
WRITE ( iw, '( A, T30, A, T37, A, T44, A, T51, A, T58, A, T71, F10.2 )' ) &
' FORCE FIELD| Jorgensen', &
ADJUSTL ( atom_names(i1)), &
ADJUSTL ( atom_names(i2)), &
ADJUSTL ( atom_names(i3)), &
ADJUSTL ( atom_names(i4)), &
torsion_parm ( i1, i2, i3, i4 ) % parm
' a0 [K]', torsion_parm ( i1, i2, i3, i4 ) % a(0)
WRITE ( iw, '( A, T58, A, T71, F10.2 )' ) ' FORCE FIELD| ', &
' a1 [K]', torsion_parm ( i1, i2, i3, i4 ) % a(1)
WRITE ( iw, '( A, T58, A, T71, F10.2 )' ) ' FORCE FIELD| ', &
' a2 [K]', torsion_parm ( i1, i2, i3, i4 ) % a(2)
WRITE ( iw, '( A, T58, A, T71, F10.2 )' ) ' FORCE FIELD| ', &
' a3 [K]', torsion_parm ( i1, i2, i3, i4 ) % a(3)
END IF
END DO
END DO
@ -828,7 +854,7 @@ SUBROUTINE input_distance_constraints ( constraint_distance, atom_names, &
DO i1 = 1, natom_types
DO i2 = i1, natom_types
IF ( constraint_distance ( i1, i2 ) /= 0.0_dbl ) THEN
WRITE ( iw, '( A, T30, A, T40, A, T71, F10.4 )' ) &
WRITE ( iw, '( A, T30, A, T37, A, T71, F10.4 )' ) &
' FORCE FIELD| Constraint: ', &
ADJUSTL ( atom_names ( i1 ) ), &
ADJUSTL ( atom_names ( i2 ) ), &

6
src/force_fit.F
View file

@ -9,7 +9,7 @@ MODULE force_fit
USE amoeba, ONLY : amoeba_evaluate
USE global_types, ONLY : globenv
USE io_parameters, ONLY : scr
! USE io_parameters, ONLY : scr
USE kinds, ONLY : dbl
USE md, ONLY : simulation_parameters_type, thermodynamic_type
USE pair_potential, ONLY : potentialparm_type, spline_nonbond_control
@ -90,11 +90,13 @@ SUBROUTINE read_fit_section ( fit_param )
TYPE ( fit_parameter_type ), INTENT ( OUT ) :: fit_param
! Locals
INTEGER :: ierror, ilen
INTEGER :: ierror, ilen, scr
CHARACTER ( len = 20 ) :: string
CHARACTER ( len = 5 ) :: label
!------------------------------------------------------------------------------
! this is just a temporary fix to get rid of io_parameters
scr = 6
!..defaults
fit_param % max_tol = 1.0E-6_dbl

19
src/global_types.F
View file

@ -2,6 +2,21 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****s* cp2k/global_types [1.0] *
!!
!! NAME
!! global_types
!!
!! FUNCTION
!!
!! AUTHOR
!! JGH
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE global_types
@ -16,6 +31,7 @@ MODULE global_types
CHARACTER ( LEN = 200 ) :: input_file_name
CHARACTER ( LEN = 200 ) :: program_name
CHARACTER ( LEN = 200 ) :: project_name
CHARACTER ( LEN = 200 ) :: default_fft_library
LOGICAL :: ionode !! this processor performs global output
INTEGER :: group !! group handle for all processors
INTEGER :: source !! processor id for global output node
@ -27,3 +43,6 @@ MODULE global_types
END TYPE global_environment_type
END MODULE global_types
!!*****
!******************************************************************************

43
src/header.F
View file

@ -2,6 +2,21 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/header [1.0] *
!!
!! NAME
!! header
!!
!! FUNCTION
!!
!! AUTHOR
!! APSI & CJM & JGH
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE header
@ -9,7 +24,11 @@ MODULE header
PRIVATE
PUBLIC :: fist_header, tbmd_header, qs_header, wave_header, faust_header
PUBLIC :: atom_header
!!*****
!******************************************************************************
CONTAINS
!******************************************************************************
@ -156,4 +175,28 @@ END SUBROUTINE faust_header
!******************************************************************************
SUBROUTINE atom_header ( iw )
IMPLICIT NONE
INTEGER :: iw
!------------------------------------------------------------------------------
WRITE ( iw, '( / )' )
WRITE ( iw, '( 14(25x,a,/) )' ) &
' **** ****** **** **** ', &
' ** ** ****** ** ** ****** ', &
' ****** ** ** ** ** ** ', &
' ** ** ** **** ** ** ', &
' ', &
' University of Zurich ', &
' 2000 ', &
' ', &
' Version 0.0 ', &
' '
END SUBROUTINE atom_header
!******************************************************************************
END MODULE header

111
src/initialize_ewalds.F
View file

@ -1,111 +0,0 @@
!-----------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
MODULE initialize_ewalds
USE dg, ONLY : dg_grid_setup
USE ewalds, ONLY : ewald_self
USE ewald_parameters_types, ONLY : ewald_parameters_type
USE io_parameters, ONLY : scr
USE kinds, ONLY : dbl
USE md, ONLY : thermodynamic_type
USE molecule_types, ONLY : particle_node_type
USE particle_types, ONLY : particle_type
USE pw_grid_types, ONLY : pw_grid_type
USE pw_grids, ONLY : pw_grid_setup
USE simulation_cell, ONLY : cell_type
PUBLIC :: initialize_ewald
PRIVATE
CONTAINS
!******************************************************************************
SUBROUTINE initialize_ewald ( part, pnode, ewald_param, box, thermo, &
ewald_grid, pme_small_grid, pme_big_grid )
IMPLICIT NONE
! Arguments
TYPE ( ewald_parameters_type ), INTENT ( OUT ) :: ewald_param
TYPE ( particle_type ), INTENT ( IN ) :: part ( : )
TYPE ( particle_node_type ), INTENT ( IN ) :: pnode ( : )
TYPE ( cell_type ), INTENT ( IN ) :: box
TYPE ( thermodynamic_type ), INTENT ( INOUT ) :: thermo
TYPE ( pw_grid_type ), INTENT ( OUT ), OPTIONAL :: ewald_grid
TYPE ( pw_grid_type ), INTENT ( OUT ), OPTIONAL :: pme_small_grid
TYPE ( pw_grid_type ), INTENT ( OUT ), OPTIONAL :: pme_big_grid
! Locals
INTEGER :: natoms, iw, iat, jat, gmax, npts ( 3 )
REAL ( dbl ) :: qi, qj, i, cutoff
!------------------------------------------------------------------------------
iw = scr
natoms = SIZE ( part )
IF ( ewald_param % ewald_type /= 'NONE' ) THEN
WRITE ( iw, '( A,T71,A )' ) ' Ewald summation is done by:', &
ewald_param % ewald_type
WRITE ( iw, '( A,T71,F10.4 )' ) ' Ewald alpha parameter [A]', &
ewald_param % alpha
SELECT CASE(ewald_param % ewald_TYPE ( 1:3))
CASE DEFAULT
WRITE ( iw, '( A,T71,I10 )' ) &
' Ewald G-space max. Miller index', ewald_param % gmax
CASE ( 'PME')
WRITE ( iw, '( A,T71,I10 )' ) &
' PME max small-grid points ', ewald_param % ns_max
WRITE ( iw, '( A,T71,F10.4 )' ) &
' PME gaussian tolerance ', ewald_param % epsilon
END SELECT
ELSE
WRITE( iw, '( A )' ) ' No Ewald summation is performed'
END IF
! fire up the reciprocal space and compute self interaction and
! term from the neutralizing background.
IF ( ewald_param % ewald_type /= 'NONE' ) THEN
CALL ewald_self ( ewald_param, thermo, pnode )
! set up EWALD_GAUSS
IF ( PRESENT ( ewald_grid)) THEN
gmax = ewald_param % gmax
IF ( gmax == 2 * ( gmax / 2 ) ) THEN
CALL stop_prg ( "initialize_ewalds", "gmax has to be odd" )
END IF
ewald_grid % bounds ( 1, : ) = -gmax / 2
ewald_grid % bounds ( 2, : ) = +gmax / 2
npts=( / gmax, gmax, gmax / )
CALL dg_find_cutoff ( npts, box, cutoff )
CALL pw_grid_setup( box, ewald_grid, cutoff)
END IF
! set up PME_GAUSS
!dbg
! IF ( PRESENT ( pme_small_grid).AND. PRESENT ( pme_big_grid))
! CALL CALL dg_setup(box, npts_s, epsilon, alpha, pme_small_grid, &
! pme_big_grid, ewald_param % ewald_type )
!dbg
END IF
END SUBROUTINE initialize_ewald
!******************************************************************************
END MODULE initialize_ewalds

62
src/initialize_extended_types.F
View file

@ -2,6 +2,21 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/initialize_extended_types [1.0] *
!!
!! NAME
!! initialize_extended_types
!!
!! FUNCTION
!!
!! AUTHOR
!! CJM
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE initialize_extended_types
@ -19,8 +34,26 @@ MODULE initialize_extended_types
PUBLIC :: initialize_extended_type
PRIVATE
!!*****
!******************************************************************************
CONTAINS
!******************************************************************************
!!****** initialize_extended_types/initialize_extended_type [1.0] *
!!
!! NAME
!! initialize_extended_types
!!
!! FUNCTION
!!
!! AUTHOR
!! CJM
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
SUBROUTINE initialize_extended_type ( box, simpar, molecule, mol_setup, nhcp, &
@ -68,6 +101,22 @@ SUBROUTINE initialize_extended_type ( box, simpar, molecule, mol_setup, nhcp, &
END SUBROUTINE initialize_extended_type
!!*****
!******************************************************************************
!!****** initialize_extended_types/get_nfree [1.0] *
!!
!! NAME
!! get_nfree
!!
!! FUNCTION
!!
!! AUTHOR
!! CJM
!!
!! MODIFICATION HISTORY
!! 29-Nov-00 (JGH) correct counting of DOF if constraints are off
!!
!! SOURCE
!******************************************************************************
SUBROUTINE get_nfree ( simpar, molecule, globenv )
@ -84,13 +133,17 @@ SUBROUTINE get_nfree ( simpar, molecule, globenv )
!------------------------------------------------------------------------------
nconst = 0
natoms = 0
DO imol = 1, SIZE ( molecule )
nconst = nconst + molecule ( imol ) % ndistconst_mol
nconst = nconst + 3 * molecule ( imol ) % ng3x3const_mol
natoms = natoms + molecule ( imol ) % natoms_mol
END DO
nconst = 0
IF ( simpar % constraint ) THEN
DO imol = 1, SIZE ( molecule )
nconst = nconst + molecule ( imol ) % ndistconst_mol
nconst = nconst + 3 * molecule ( imol ) % ng3x3const_mol
END DO
END IF
#if defined(__parallel)
CALL mp_sum(nconst,globenv % group)
@ -115,6 +168,9 @@ SUBROUTINE get_nfree ( simpar, molecule, globenv )
END SUBROUTINE get_nfree
!!*****
!******************************************************************************
END MODULE initialize_extended_types
!******************************************************************************

36
src/initialize_molecule_types.F
View file

@ -2,6 +2,21 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/initialize_molecule_types [1.0] *
!!
!! NAME
!! initialize_molecule_types
!!
!! FUNCTION
!!
!! AUTHOR
!! CJM
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE initialize_molecule_types
@ -20,8 +35,26 @@ MODULE initialize_molecule_types
PUBLIC :: initialize_molecule_type
PRIVATE
!!*****
!******************************************************************************
CONTAINS
!******************************************************************************
!!****** initialize_molecule_types/initialize_molecule_type [1.0] *
!!
!! NAME
!! initialize_molecule_types
!!
!! FUNCTION
!!
!! AUTHOR
!! CJM
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
SUBROUTINE initialize_molecule_type ( mol_setup, intra_param, &
@ -135,6 +168,9 @@ SUBROUTINE initialize_molecule_type ( mol_setup, intra_param, &
END SUBROUTINE initialize_molecule_type
!!*****
!******************************************************************************
END MODULE initialize_molecule_types
!******************************************************************************

34
src/initialize_particle_types.F
View file

@ -2,6 +2,20 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/initialize_particle_types [1.0] *
!!
!! NAME
!! initialize_particle_types
!!
!! FUNCTION
!!
!! AUTHOR
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE initialize_particle_types
@ -18,8 +32,25 @@ MODULE initialize_particle_types
PUBLIC :: initialize_particle_type
PRIVATE
!!*****
!******************************************************************************
CONTAINS
!******************************************************************************
!!****** initialize_particle_types/initialize_particle_type [1.0] *
!!
!! NAME
!! initialize_particle_type
!!
!! FUNCTION
!!
!! AUTHOR
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
SUBROUTINE initialize_particle_type ( atom_names, simpar, mol_setup, &
@ -76,6 +107,9 @@ SUBROUTINE initialize_particle_type ( atom_names, simpar, mol_setup, &
END SUBROUTINE initialize_particle_type
!!*****
!******************************************************************************
END MODULE initialize_particle_types
!******************************************************************************

18
src/input_types.F
View file

@ -2,6 +2,21 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****s* cp2k/input_types [1.0] *
!!
!! NAME
!! input_types
!!
!! FUNCTION
!!
!! AUTHOR
!! CJM & JGH
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE input_types
@ -21,3 +36,6 @@ MODULE input_types
PUBLIC :: setup_parameters_type
END MODULE input_types
!!*****
!******************************************************************************

105
src/integrator.F
View file

@ -2,19 +2,33 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/integrator [1.0] *
!!
!! NAME
!! integrator
!!
!! FUNCTION
!! Provides integrator routines (velocity verlet) for all the
!! ensemble types
!!
!! AUTHOR
!! CJM
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE integrator
! USE constraint, ONLY : shake_control, rattle_control, getold, &
! pv_constraint, initialize_roll, rattle_roll_control, &
! shake_roll_control
USE constraint
USE constraint, ONLY: shake_control, rattle_control, pv_constraint, getold, &
roll, initialize_roll, shake_roll_control, rattle_roll_control
USE global_types, ONLY: global_environment_type
USE eigenvalueproblems, ONLY : diagonalise
USE ewald_parameters_types, ONLY : ewald_parameters_type
USE force_control, ONLY : force
USE kinds, ONLY : dbl
!*apsi USE lapack, ONLY : rs
USE linklist_control, ONLY : get_ll_parm
USE mathconstants, ONLY : pi, zero
USE md, ONLY : simulation_parameters_type, thermodynamic_type, &
@ -67,10 +81,26 @@ MODULE integrator
TYPE ( global_environment_type ) :: intenv
!!*****
!******************************************************************************
CONTAINS
!******************************************************************************
!!****** integrator/velocity_verlet [1.0] *
!!
!! NAME
!! velocity_verlet
!!
!! FUNCTION
!!
!! AUTHOR
!! CJM
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
SUBROUTINE velocity_verlet ( itimes, constant, simpar, inter, &
@ -113,10 +143,29 @@ SUBROUTINE velocity_verlet ( itimes, constant, simpar, inter, &
END SUBROUTINE velocity_verlet
!!*****
!******************************************************************************
!!****** integrator/set_integrator [1.0] *
!!
!! NAME
!! set_integrator
!!
!! FUNCTION
!! setup of parameters for integrator
!!
!! AUTHOR
!! JGH
!!
!! MODIFICATION HISTORY
!! none
!!
!! NOTES
!! Call this subroutine before the first call to energy or velocity_verlet
!! or if you want to change ionode and/or output files
!!
!! SOURCE
!******************************************************************************
! call this subroutine before the first call to energy or velocity_verlet
! or if you want to change ionode and/or output files
SUBROUTINE set_integrator ( globenv, mdio )
IMPLICIT NONE
@ -150,6 +199,22 @@ SUBROUTINE set_integrator ( globenv, mdio )
END SUBROUTINE set_integrator
!!*****
!******************************************************************************
!!****** integrator/set_energy_parm [1.0] *
!!
!! NAME
!! set_energy_parm
!!
!! FUNCTION
!!
!! AUTHOR
!! CJM
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
SUBROUTINE set_energy_parm ( pconv, tempconv, l, vol, e, pres, temp, angl )
@ -173,9 +238,26 @@ SUBROUTINE set_energy_parm ( pconv, tempconv, l, vol, e, pres, temp, angl )
IF ( PRESENT ( angl)) angl_label = angl
END SUBROUTINE set_energy_parm
!!*****
!******************************************************************************
! this routine computes the conserved quantity, temperature
! and things like that and prints them out
!!****** integrator/energy [1.0] *
!!
!! NAME
!! energy
!!
!! FUNCTION
!! This routine computes the conserved quantity, temperature
!! and things like that and prints them out
!!
!! AUTHOR
!! CJM
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
SUBROUTINE energy ( itimes, constant, simpar, struc, thermo, nhcp )
@ -418,7 +500,7 @@ IF (ionode .AND. mod(itimes,ivel)==0) THEN
adjustl(temp_label), '= ', temp, averages%avetemp
WRITE ( iw, '( A,A,T31,A,T49,2(1X,E15.7) )' ) ' PRESSURE', &
adjustl(pres_label), '= ', pv_scalar, averages%avepress
CALL get_ll_parm(lup,aup)
CALL get_ll_parm(struc%ll_data,lup,aup)
IF (simpar%ensemble(1:3)=='NPT') THEN
WRITE ( iw, '( A,A,T31,A,T49,2(1X,E15.7) )' ) ' VOLUME', &
adjustl(vol_label), '= ', struc%box%deth, averages%avevol
@ -446,6 +528,7 @@ IF (ionode .AND. mod(itimes,ivel)==0) THEN
END SUBROUTINE energy
!!*****
!******************************************************************************
! nve integrator for particle positions & momenta
@ -1114,3 +1197,5 @@ END SUBROUTINE rattle_init
!******************************************************************************
END MODULE integrator
!******************************************************************************

55
src/io_parameters.F
View file

@ -1,55 +0,0 @@
!-----------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
MODULE io_parameters
IMPLICIT NONE
! Defaults for some global variables
! *****************************************************************************
! scr: Logical unit number of the screen
! crd: Logical unit number of the coordinate file
! vel: Logical unit number of the velocity file
! ptn: Logical unit number of the pressure tensor file
! ene: Logical unit number of the energy file
! tem: Logical unit number of the temperature file
! dmp: Logical unit number of the dump file
! icrd : frequency of writing coordinate file
! ivel : frequency of writing velocity file
! iptens : frequency of writing pressure tensor file
! itemp : frequency of writing temperature file
! idump : frequency of writing dump file
! iscreen : frequency of writing screen file
! *****************************************************************************
PRIVATE
PUBLIC :: crd_file_name, vel_file_name, ptens_file_name, ener_file_name, &
temp_file_name, dump_file_name
PUBLIC :: crd, vel, ptn, ene, tem, scr, dmp
PUBLIC :: icrd, ivel, iptens, iener, itemp, iscreen, idump
CHARACTER ( LEN = 50 ) :: crd_file_name, &
vel_file_name, &
ptens_file_name,&
ener_file_name, &
temp_file_name, &
dump_file_name
INTEGER, PARAMETER :: crd = 1, &
vel = 2, &
ptn = 3, &
ene = 4, &
tem = 5, &
scr = 6, &
dmp = 7
INTEGER :: icrd, ivel, iptens, &
iener, itemp, idump, &
iscreen
END MODULE io_parameters

24
src/k290.F
View file

@ -2,8 +2,26 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!------------------------------------------------------------------------------!
! invmat uses lapack routines --> should be handled by a more general f95
! interface
!!****** cp2k/k290 [1.0] *
!!
!! NAME
!! k290
!!
!! FUNCTION
!! K-points and symmetry routines
!!
!! AUTHOR
!! JGH
!!
!! MODIFICATION HISTORY
!! none
!!
!! NOTES
!! invmat uses lapack routines --> should be handled by a more general f95
!! interface
!!
!! SOURCE
!******************************************************************************
!
MODULE k290
!
@ -27,6 +45,8 @@
REAL (dbl) :: a1(3), a2(3), a3(3), alat
REAL (dbl) :: strain(3,3) = 0._dbl, wvk0(3) = 0._dbl
REAL (dbl) :: delta = 1.E-6_dbl
!!*****
!------------------------------------------------------------------------------!
!
CONTAINS

58
src/kinds.F
View file

@ -2,14 +2,32 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!! Data type definitions; tested on:
!! - IBM AIX xlf90
!! - SGI IRIX f90
!! - CRAY T3E f90
!! - DEC ALPHA f90
!! - NAG_F90
!! - SUN
!! - HITACHI
!!****** cp2k/kinds [1.0] *
!!
!! NAME
!! kinds
!!
!! FUNCTION
!! Defines the basic variable types
!!
!! AUTHOR
!! Matthias Krack
!!
!! MODIFICATION HISTORY
!! Adapted for CP2K by JGH
!!
!! NOTES
!! Data type definitions; tested on:
!! - IBM AIX xlf90
!! - SGI IRIX f90
!! - CRAY T3E f90
!! - DEC ALPHA f90
!! - NAG_F90
!! - SUN
!! - HITACHI
!!
!! SOURCE
!******************************************************************************
MODULE kinds
@ -26,12 +44,29 @@ MODULE kinds
LOGICAL, PARAMETER :: YES = .TRUE., NO = .FALSE.
LOGICAL, PARAMETER :: TRUE = .TRUE., FALSE = .FALSE.
!!*****
!******************************************************************************
CONTAINS
!******************************************************************************
!! Print informations about the used data types.
!!****** kinds/print_kind_info [1.0] *
!!
!! NAME
!! print_kind_info
!!
!! FUNCTION
!! Print informations about the used data types.
!!
!! AUTHOR
!! Matthias Krack
!!
!! MODIFICATION HISTORY
!! Adapted by JGH for Cp2k
!!
!! SOURCE
!******************************************************************************
SUBROUTINE print_kind_info ( iw )
@ -72,6 +107,9 @@ SUBROUTINE print_kind_info ( iw )
END SUBROUTINE print_kind_info
!!*****
!******************************************************************************
END MODULE kinds
!******************************************************************************

37
src/kpoint_initialization.F
View file

@ -2,6 +2,21 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!------------------------------------------------------------------------------!
!!****** cp2k/kpoint_initialization [1.0] *
!!
!! NAME
!! kpoint_initialization
!!
!! FUNCTION
!!
!! AUTHOR
!! JGH
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
!
MODULE kpoint_initialization
!
@ -18,11 +33,29 @@
PRIVATE
!
PUBLIC :: initialize_kpoints
!!*****
!------------------------------------------------------------------------------!
!
CONTAINS
!
!------------------------------------------------------------------------------!
!!****** kpoint_initialization/initialize_kpoints [1.0] *
!!
!! NAME
!! initialize_kpoints
!!
!! FUNCTION
!!
!! AUTHOR
!! JGH
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
SUBROUTINE initialize_kpoints(inpar,kp,symm,hmat,part)
IMPLICIT NONE
@ -59,6 +92,10 @@
END IF
END SUBROUTINE initialize_kpoints
!!*****
!------------------------------------------------------------------------------!
END MODULE kpoint_initialization
!------------------------------------------------------------------------------!

15
src/lapack.F
View file

@ -2,6 +2,20 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/lapack [1.0] *
!!
!! NAME
!! lapack
!!
!! FUNCTION
!!
!! AUTHOR
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
#if ! defined ( __YMP ) && ! defined ( __T3E )
@ -120,3 +134,4 @@ MODULE lapack
END MODULE lapack
#endif
!!*****

2
src/lib/Makefile
View file

@ -20,6 +20,6 @@ clean:
rm -f *.o *.mod
realclean:
rm -f *.a *.o *.mod *~ F*.f *.lst
rm -f $(LIB) *.o *.mod *~ F*.f *.lst
####################

56
src/lib/fftpre.F
View file

@ -24,7 +24,7 @@ SUBROUTINE fftpre ( mm, nfft, m, nn, n, zin, zout, &
r, r1, r2, r3, r4, r5, r6, r7, r8, cr2, cr3, cr4, cr5, &
ci2, ci3, ci4, ci5, ur1, ur2, ur3, ui1, ui2, ui3, &
vr1, vr2, vr3, vi1, vi2, vi3, cm, cp, dm, dp, &
am, ap, bm, bp,s25, s34, r34, r25, sin2, sin4
am, ap, bm, bp,s25, s34, r34, r25, sin2, sin4, bbs
REAL ( dbl ), PARAMETER :: rt2i = 0.7071067811865475_dbl ! sqrt(0.5)
REAL ( dbl ), PARAMETER :: bb = 0.8660254037844387_dbl ! sqrt(3)/2
REAL ( dbl ), PARAMETER :: cos2 = 0.3090169943749474_dbl ! cos(2*pi/5)
@ -510,6 +510,7 @@ SUBROUTINE fftpre ( mm, nfft, m, nn, n, zin, zout, &
END DO
END IF
ELSE IF ( now == 3 ) THEN
bbs = bb * REAL ( isign, dbl )
ia = 1
nin1 = ia - after
nout1 = ia - atn
@ -533,8 +534,8 @@ SUBROUTINE fftpre ( mm, nfft, m, nn, n, zin, zout, &
zout ( 2, j, nout1 ) = s + s1
r1 = r1 - 0.5_dbl * r
s1 = s1 - 0.5_dbl * s
r2 = bb * ( r2 - r3 )
s2 = bb * ( s2 - s3 )
r2 = bbs * ( r2 - r3 )
s2 = bbs * ( s2 - s3 )
zout ( 1, j, nout2 ) = r1 - s2
zout ( 2, j, nout2 ) = s1 + r2
zout ( 1, j, nout3 ) = r1 + s2
@ -566,9 +567,9 @@ SUBROUTINE fftpre ( mm, nfft, m, nn, n, zin, zout, &
zout ( 1, j, nout1 ) = r + r1
zout ( 2, j, nout1 ) = s + s1
r1 = r1 - 0.5_dbl * r
s1 = s1 - 0.5_dbl * s
r2 = bb*(r2-r3)
s2 = bb*(s2-s3)
s1 = s1 - 0.5_dbl * s
r2 = bbs*(r2-r3)
s2 = bbs*(s2-s3)
zout ( 1, j, nout2 ) = r1 - s2
zout ( 2, j, nout2 ) = s1 + r2
zout ( 1, j, nout3 ) = r1 + s2
@ -597,9 +598,9 @@ SUBROUTINE fftpre ( mm, nfft, m, nn, n, zin, zout, &
zout ( 1, j, nout1 ) = r + r1
zout ( 2, j, nout1 ) = s + s1
r1 = r1 - 0.5_dbl * r
s1 = s1 - 0.5_dbl * s
r2 = bb * ( r2 - r3 )
s2 = bb * ( s2 - s3 )
s1 = s1 - 0.5_dbl * s
r2 = bbs * ( r2 - r3 )
s2 = bbs * ( s2 - s3 )
zout ( 1, j, nout2 ) = r1 - s2
zout ( 2, j, nout2 ) = s1 + r2
zout ( 1, j, nout3 ) = r1 + s2
@ -632,9 +633,9 @@ SUBROUTINE fftpre ( mm, nfft, m, nn, n, zin, zout, &
zout ( 1, j, nout1 ) = r + r1
zout ( 2, j, nout1 ) = s + s1
r1 = r1 - 0.5_dbl * r
s1 = s1 - 0.5_dbl * s
r2 = bb * ( r2 - r3 )
s2 = bb * ( s2 - s3 )
s1 = s1 - 0.5_dbl * s
r2 = bbs * ( r2 - r3 )
s2 = bbs * ( s2 - s3 )
zout ( 1, j, nout2 ) = r1 - s2
zout ( 2, j, nout2 ) = s1 + r2
zout ( 1, j, nout3 ) = r1 + s2
@ -665,9 +666,9 @@ SUBROUTINE fftpre ( mm, nfft, m, nn, n, zin, zout, &
zout ( 1, j, nout1 ) = r + r1
zout ( 2, j, nout1 ) = s + s1
r1 = r1 - 0.5_dbl * r
s1 = s1 - 0.5_dbl * s
r2 = bb * ( r2 - r3 )
s2 = bb * ( s2 - s3 )
s1 = s1 - 0.5_dbl * s
r2 = bbs * ( r2 - r3 )
s2 = bbs * ( s2 - s3 )
zout ( 1, j, nout2 ) = r1 - s2
zout ( 2, j, nout2 ) = s1 + r2
zout ( 1, j, nout3 ) = r1 + s2
@ -709,8 +710,8 @@ SUBROUTINE fftpre ( mm, nfft, m, nn, n, zin, zout, &
zout ( 2, j, nout1 ) = s + s1
r1 = r1 - 0.5_dbl * r
s1 = s1 - 0.5_dbl * s
r2 = bb * ( r2 - r3 )
s2 = bb * ( s2 - s3 )
r2 = bbs * ( r2 - r3 )
s2 = bbs * ( s2 - s3 )
zout ( 1, j, nout2 ) = r1 - s2
zout ( 2, j, nout2 ) = s1 + r2
zout ( 1, j, nout3 ) = r1 + s2
@ -786,7 +787,7 @@ SUBROUTINE fftpre ( mm, nfft, m, nn, n, zin, zout, &
nin2 = nin1 + atb
nin3 = nin2 + atb
nin4 = nin3 + atb
nin5 = nin4 + atb
nin5 = nin4 + atb
nout1 = nout1 + atn
nout2 = nout1 + after
nout3 = nout2 + after
@ -843,7 +844,7 @@ SUBROUTINE fftpre ( mm, nfft, m, nn, n, zin, zout, &
nin2 = nin1 + atb
nin3 = nin2 + atb
nin4 = nin3 + atb
nin5 = nin4 + atb
nin5 = nin4 + atb
nout1 = nout1 + atn
nout2 = nout1 + after
nout3 = nout2 + after
@ -971,6 +972,7 @@ SUBROUTINE fftpre ( mm, nfft, m, nn, n, zin, zout, &
END IF
END DO
ELSE IF ( now == 6 ) THEN
bbs = bb * REAL ( isign, dbl )
ia = 1
nin1 = ia - after
nout1 = ia - atn
@ -1002,10 +1004,10 @@ SUBROUTINE fftpre ( mm, nfft, m, nn, n, zin, zout, &
s1 = s1 - 0.5_dbl * s
r = r2 - r3
s = s2 - s3
ur2 = r1 - s * bb
ui2 = s1 + r * bb
ur3 = r1 + s * bb
ui3 = s1 - r * bb
ur2 = r1 - s * bbs
ui2 = s1 + r * bbs
ur3 = r1 + s * bbs
ui3 = s1 - r * bbs
r2 = zin ( 1, nin6, j )
s2 = zin ( 2, nin6, j )
@ -1021,10 +1023,10 @@ SUBROUTINE fftpre ( mm, nfft, m, nn, n, zin, zout, &
s1 = s1 - 0.5_dbl * s
r = r2 - r3
s = s2 - s3
vr2 = r1 - s * bb
vi2 = s1 + r * bb
vr3 = r1 + s * bb
vi3 = s1 - r * bb
vr2 = r1 - s * bbs
vi2 = s1 + r * bbs
vr3 = r1 + s * bbs
vi3 = s1 - r * bbs
zout ( 1, j, nout1 ) = ur1 + vr1
zout ( 2, j, nout1 ) = ui1 + vi1

56
src/lib/fftrot.F
View file

@ -24,7 +24,7 @@ SUBROUTINE fftrot ( mm, nfft, m, nn, n, zin, zout, &
r, r1, r2, r3, r4, r5, r6, r7, r8, cr2, cr3, cr4, cr5, &
ci2, ci3, ci4, ci5, ur1, ur2, ur3, ui1, ui2, ui3, &
vr1, vr2, vr3, vi1, vi2, vi3, cm, cp, dm, dp, &
am, ap, bm, bp,s25, s34, r34, r25, sin2, sin4
am, ap, bm, bp,s25, s34, r34, r25, sin2, sin4, bbs
REAL ( dbl ), PARAMETER :: rt2i = 0.7071067811865475_dbl ! sqrt(0.5)
REAL ( dbl ), PARAMETER :: bb = 0.8660254037844387_dbl ! sqrt(3)/2
REAL ( dbl ), PARAMETER :: cos2 = 0.3090169943749474_dbl ! cos(2*pi/5)
@ -510,6 +510,7 @@ SUBROUTINE fftrot ( mm, nfft, m, nn, n, zin, zout, &
END DO
END IF
ELSE IF ( now == 3 ) THEN
bbs = bb * REAL ( isign, dbl )
ia = 1
nin1 = ia - after
nout1 = ia - atn
@ -533,8 +534,8 @@ SUBROUTINE fftrot ( mm, nfft, m, nn, n, zin, zout, &
zout ( 2, nout1, j ) = s + s1
r1 = r1 - 0.5_dbl * r
s1 = s1 - 0.5_dbl * s
r2 = bb * ( r2 - r3 )
s2 = bb * ( s2 - s3 )
r2 = bbs * ( r2 - r3 )
s2 = bbs * ( s2 - s3 )
zout ( 1, nout2, j ) = r1 - s2
zout ( 2, nout2, j ) = s1 + r2
zout ( 1, nout3, j ) = r1 + s2
@ -566,9 +567,9 @@ SUBROUTINE fftrot ( mm, nfft, m, nn, n, zin, zout, &
zout ( 1, nout1, j ) = r + r1
zout ( 2, nout1, j ) = s + s1
r1 = r1 - 0.5_dbl * r
s1 = s1 - 0.5_dbl * s
r2 = bb*(r2-r3)
s2 = bb*(s2-s3)
s1 = s1 - 0.5_dbl * s
r2 = bbs*(r2-r3)
s2 = bbs*(s2-s3)
zout ( 1, nout2, j ) = r1 - s2
zout ( 2, nout2, j ) = s1 + r2
zout ( 1, nout3, j ) = r1 + s2
@ -597,9 +598,9 @@ SUBROUTINE fftrot ( mm, nfft, m, nn, n, zin, zout, &
zout ( 1, nout1, j ) = r + r1
zout ( 2, nout1, j ) = s + s1
r1 = r1 - 0.5_dbl * r
s1 = s1 - 0.5_dbl * s
r2 = bb * ( r2 - r3 )
s2 = bb * ( s2 - s3 )
s1 = s1 - 0.5_dbl * s
r2 = bbs * ( r2 - r3 )
s2 = bbs * ( s2 - s3 )
zout ( 1, nout2, j ) = r1 - s2
zout ( 2, nout2, j ) = s1 + r2
zout ( 1, nout3, j ) = r1 + s2
@ -632,9 +633,9 @@ SUBROUTINE fftrot ( mm, nfft, m, nn, n, zin, zout, &
zout ( 1, nout1, j ) = r + r1
zout ( 2, nout1, j ) = s + s1
r1 = r1 - 0.5_dbl * r
s1 = s1 - 0.5_dbl * s
r2 = bb * ( r2 - r3 )
s2 = bb * ( s2 - s3 )
s1 = s1 - 0.5_dbl * s
r2 = bbs * ( r2 - r3 )
s2 = bbs * ( s2 - s3 )
zout ( 1, nout2, j ) = r1 - s2
zout ( 2, nout2, j ) = s1 + r2
zout ( 1, nout3, j ) = r1 + s2
@ -665,9 +666,9 @@ SUBROUTINE fftrot ( mm, nfft, m, nn, n, zin, zout, &
zout ( 1, nout1, j ) = r + r1
zout ( 2, nout1, j ) = s + s1
r1 = r1 - 0.5_dbl * r
s1 = s1 - 0.5_dbl * s
r2 = bb * ( r2 - r3 )
s2 = bb * ( s2 - s3 )
s1 = s1 - 0.5_dbl * s
r2 = bbs * ( r2 - r3 )
s2 = bbs * ( s2 - s3 )
zout ( 1, nout2, j ) = r1 - s2
zout ( 2, nout2, j ) = s1 + r2
zout ( 1, nout3, j ) = r1 + s2
@ -709,8 +710,8 @@ SUBROUTINE fftrot ( mm, nfft, m, nn, n, zin, zout, &
zout ( 2, nout1, j ) = s + s1
r1 = r1 - 0.5_dbl * r
s1 = s1 - 0.5_dbl * s
r2 = bb * ( r2 - r3 )
s2 = bb * ( s2 - s3 )
r2 = bbs * ( r2 - r3 )
s2 = bbs * ( s2 - s3 )
zout ( 1, nout2, j ) = r1 - s2
zout ( 2, nout2, j ) = s1 + r2
zout ( 1, nout3, j ) = r1 + s2
@ -786,7 +787,7 @@ SUBROUTINE fftrot ( mm, nfft, m, nn, n, zin, zout, &
nin2 = nin1 + atb
nin3 = nin2 + atb
nin4 = nin3 + atb
nin5 = nin4 + atb
nin5 = nin4 + atb
nout1 = nout1 + atn
nout2 = nout1 + after
nout3 = nout2 + after
@ -843,7 +844,7 @@ SUBROUTINE fftrot ( mm, nfft, m, nn, n, zin, zout, &
nin2 = nin1 + atb
nin3 = nin2 + atb
nin4 = nin3 + atb
nin5 = nin4 + atb
nin5 = nin4 + atb
nout1 = nout1 + atn
nout2 = nout1 + after
nout3 = nout2 + after
@ -971,6 +972,7 @@ SUBROUTINE fftrot ( mm, nfft, m, nn, n, zin, zout, &
END IF
END DO
ELSE IF ( now == 6 ) THEN
bbs = bb * REAL ( isign, dbl )
ia = 1
nin1 = ia - after
nout1 = ia - atn
@ -1002,10 +1004,10 @@ SUBROUTINE fftrot ( mm, nfft, m, nn, n, zin, zout, &
s1 = s1 - 0.5_dbl * s
r = r2 - r3
s = s2 - s3
ur2 = r1 - s * bb
ui2 = s1 + r * bb
ur3 = r1 + s * bb
ui3 = s1 - r * bb
ur2 = r1 - s * bbs
ui2 = s1 + r * bbs
ur3 = r1 + s * bbs
ui3 = s1 - r * bbs
r2 = zin ( 1, j, nin6 )
s2 = zin ( 2, j, nin6 )
@ -1021,10 +1023,10 @@ SUBROUTINE fftrot ( mm, nfft, m, nn, n, zin, zout, &
s1 = s1 - 0.5_dbl * s
r = r2 - r3
s = s2 - s3
vr2 = r1 - s * bb
vi2 = s1 + r * bb
vr3 = r1 + s * bb
vi3 = s1 - r * bb
vr2 = r1 - s * bbs
vi2 = s1 + r * bbs
vr3 = r1 + s * bbs
vi3 = s1 - r * bbs
zout ( 1, nout1, j ) = ur1 + vr1
zout ( 2, nout1, j ) = ui1 + vi1

56
src/lib/fftstp.F
View file

@ -24,7 +24,7 @@ SUBROUTINE fftstp ( mm, nfft, m, nn, n, zin, zout, &
r, r1, r2, r3, r4, r5, r6, r7, r8, cr2, cr3, cr4, cr5, &
ci2, ci3, ci4, ci5, ur1, ur2, ur3, ui1, ui2, ui3, &
vr1, vr2, vr3, vi1, vi2, vi3, cm, cp, dm, dp, &
am, ap, bm, bp,s25, s34, r34, r25, sin2, sin4
am, ap, bm, bp,s25, s34, r34, r25, sin2, sin4, bbs
REAL ( dbl ), PARAMETER :: rt2i = 0.7071067811865475_dbl ! sqrt(0.5)
REAL ( dbl ), PARAMETER :: bb = 0.8660254037844387_dbl ! sqrt(3)/2
REAL ( dbl ), PARAMETER :: cos2 = 0.3090169943749474_dbl ! cos(2*pi/5)
@ -510,6 +510,7 @@ SUBROUTINE fftstp ( mm, nfft, m, nn, n, zin, zout, &
END DO
END IF
ELSE IF ( now == 3 ) THEN
bbs = bb * REAL ( isign, dbl )
ia = 1
nin1 = ia - after
nout1 = ia - atn
@ -533,8 +534,8 @@ SUBROUTINE fftstp ( mm, nfft, m, nn, n, zin, zout, &
zout ( 2, j, nout1 ) = s + s1
r1 = r1 - 0.5_dbl * r
s1 = s1 - 0.5_dbl * s
r2 = bb * ( r2 - r3 )
s2 = bb * ( s2 - s3 )
r2 = bbs * ( r2 - r3 )
s2 = bbs * ( s2 - s3 )
zout ( 1, j, nout2 ) = r1 - s2
zout ( 2, j, nout2 ) = s1 + r2
zout ( 1, j, nout3 ) = r1 + s2
@ -566,9 +567,9 @@ SUBROUTINE fftstp ( mm, nfft, m, nn, n, zin, zout, &
zout ( 1, j, nout1 ) = r + r1
zout ( 2, j, nout1 ) = s + s1
r1 = r1 - 0.5_dbl * r
s1 = s1 - 0.5_dbl * s
r2 = bb*(r2-r3)
s2 = bb*(s2-s3)
s1 = s1 - 0.5_dbl * s
r2 = bbs*(r2-r3)
s2 = bbs*(s2-s3)
zout ( 1, j, nout2 ) = r1 - s2
zout ( 2, j, nout2 ) = s1 + r2
zout ( 1, j, nout3 ) = r1 + s2
@ -597,9 +598,9 @@ SUBROUTINE fftstp ( mm, nfft, m, nn, n, zin, zout, &
zout ( 1, j, nout1 ) = r + r1
zout ( 2, j, nout1 ) = s + s1
r1 = r1 - 0.5_dbl * r
s1 = s1 - 0.5_dbl * s
r2 = bb * ( r2 - r3 )
s2 = bb * ( s2 - s3 )
s1 = s1 - 0.5_dbl * s
r2 = bbs * ( r2 - r3 )
s2 = bbs * ( s2 - s3 )
zout ( 1, j, nout2 ) = r1 - s2
zout ( 2, j, nout2 ) = s1 + r2
zout ( 1, j, nout3 ) = r1 + s2
@ -632,9 +633,9 @@ SUBROUTINE fftstp ( mm, nfft, m, nn, n, zin, zout, &
zout ( 1, j, nout1 ) = r + r1
zout ( 2, j, nout1 ) = s + s1
r1 = r1 - 0.5_dbl * r
s1 = s1 - 0.5_dbl * s
r2 = bb * ( r2 - r3 )
s2 = bb * ( s2 - s3 )
s1 = s1 - 0.5_dbl * s
r2 = bbs * ( r2 - r3 )
s2 = bbs * ( s2 - s3 )
zout ( 1, j, nout2 ) = r1 - s2
zout ( 2, j, nout2 ) = s1 + r2
zout ( 1, j, nout3 ) = r1 + s2
@ -665,9 +666,9 @@ SUBROUTINE fftstp ( mm, nfft, m, nn, n, zin, zout, &
zout ( 1, j, nout1 ) = r + r1
zout ( 2, j, nout1 ) = s + s1
r1 = r1 - 0.5_dbl * r
s1 = s1 - 0.5_dbl * s
r2 = bb * ( r2 - r3 )
s2 = bb * ( s2 - s3 )
s1 = s1 - 0.5_dbl * s
r2 = bbs * ( r2 - r3 )
s2 = bbs * ( s2 - s3 )
zout ( 1, j, nout2 ) = r1 - s2
zout ( 2, j, nout2 ) = s1 + r2
zout ( 1, j, nout3 ) = r1 + s2
@ -709,8 +710,8 @@ SUBROUTINE fftstp ( mm, nfft, m, nn, n, zin, zout, &
zout ( 2, j, nout1 ) = s + s1
r1 = r1 - 0.5_dbl * r
s1 = s1 - 0.5_dbl * s
r2 = bb * ( r2 - r3 )
s2 = bb * ( s2 - s3 )
r2 = bbs * ( r2 - r3 )
s2 = bbs * ( s2 - s3 )
zout ( 1, j, nout2 ) = r1 - s2
zout ( 2, j, nout2 ) = s1 + r2
zout ( 1, j, nout3 ) = r1 + s2
@ -786,7 +787,7 @@ SUBROUTINE fftstp ( mm, nfft, m, nn, n, zin, zout, &
nin2 = nin1 + atb
nin3 = nin2 + atb
nin4 = nin3 + atb
nin5 = nin4 + atb
nin5 = nin4 + atb
nout1 = nout1 + atn
nout2 = nout1 + after
nout3 = nout2 + after
@ -843,7 +844,7 @@ SUBROUTINE fftstp ( mm, nfft, m, nn, n, zin, zout, &
nin2 = nin1 + atb
nin3 = nin2 + atb
nin4 = nin3 + atb
nin5 = nin4 + atb
nin5 = nin4 + atb
nout1 = nout1 + atn
nout2 = nout1 + after
nout3 = nout2 + after
@ -971,6 +972,7 @@ SUBROUTINE fftstp ( mm, nfft, m, nn, n, zin, zout, &
END IF
END DO
ELSE IF ( now == 6 ) THEN
bbs = bb * REAL ( isign, dbl )
ia = 1
nin1 = ia - after
nout1 = ia - atn
@ -1002,10 +1004,10 @@ SUBROUTINE fftstp ( mm, nfft, m, nn, n, zin, zout, &
s1 = s1 - 0.5_dbl * s
r = r2 - r3
s = s2 - s3
ur2 = r1 - s * bb
ui2 = s1 + r * bb
ur3 = r1 + s * bb
ui3 = s1 - r * bb
ur2 = r1 - s * bbs
ui2 = s1 + r * bbs
ur3 = r1 + s * bbs
ui3 = s1 - r * bbs
r2 = zin ( 1, j, nin6 )
s2 = zin ( 2, j, nin6 )
@ -1021,10 +1023,10 @@ SUBROUTINE fftstp ( mm, nfft, m, nn, n, zin, zout, &
s1 = s1 - 0.5_dbl * s
r = r2 - r3
s = s2 - s3
vr2 = r1 - s * bb
vi2 = s1 + r * bb
vr3 = r1 + s * bb
vi3 = s1 - r * bb
vr2 = r1 - s * bbs
vi2 = s1 + r * bbs
vr3 = r1 + s * bbs
vi3 = s1 - r * bbs
zout ( 1, j, nout1 ) = ur1 + vr1
zout ( 2, j, nout1 ) = ui1 + vi1

214
src/library_tests.F
View file

@ -2,6 +2,22 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/library_tests [1.0] *
!!
!! NAME
!! library_tests
!!
!! FUNCTION
!! Performance tests for basic tasks like matrix multiplies, copy, fft.
!!
!! AUTHOR
!! JGH 6-NOV-2000
!!
!! MODIFICATION HISTORY
!! 30-Nov-2000 (JGH) added input
!!
!! SOURCE
!******************************************************************************
MODULE library_tests
@ -11,6 +27,9 @@ MODULE library_tests
USE global_types, ONLY : global_environment_type
USE machine, ONLY : m_cputime
USE stop_program, ONLY : stop_memory
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
IMPLICIT NONE
@ -18,11 +37,29 @@ MODULE library_tests
PUBLIC :: lib_test
INTEGER :: runtest ( 100 )
REAL ( dbl ) :: max_memory
!!*****
!******************************************************************************
CONTAINS
!******************************************************************************
!!****** library_tests/lib_test [1.0] *
!!
!! NAME
!! lib_test
!!
!! FUNCTION
!! Master routine for tests
!!
!! AUTHOR
!! JGH 6-NOV-2000
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
SUBROUTINE lib_test ( globenv )
@ -38,18 +75,115 @@ SUBROUTINE lib_test ( globenv )
WRITE ( iw, '(T2,79("*"))' )
END IF
!
! CALL test_input ( globenv )
runtest = 0
runtest ( 3 ) = 1
CALL test_input ( globenv )
!
IF ( runtest ( 1 ) == 1 ) CALL copy_test ( globenv )
IF ( runtest ( 1 ) /= 0 ) CALL copy_test ( globenv )
!
IF ( runtest ( 2 ) == 1 ) CALL matmul_test ( globenv )
IF ( runtest ( 2 ) /= 0 ) CALL matmul_test ( globenv )
!
IF ( runtest ( 3 ) == 1 ) CALL fft_test ( globenv )
IF ( runtest ( 3 ) /= 0 ) CALL fft_test ( globenv )
END SUBROUTINE lib_test
!!*****
!******************************************************************************
!!****** library_tests/test_input [1.0] *
!!
!! NAME
!! test_input
!!
!! FUNCTION
!! Reads input section &TEST ... &END
!!
!! AUTHOR
!! JGH 30-NOV-2000
!!
!! MODIFICATION HISTORY
!! none
!!
!! NOTES
!!I---------------------------------------------------------------------------I
!!I SECTION: &TEST ... &END I
!!I I
!!I MEMORY max_memory I
!!I COPY n I
!!I MATMUL n I
!!I FFT n I
!!I I
!!I---------------------------------------------------------------------------I
!!
!! SOURCE
!******************************************************************************
SUBROUTINE test_input ( globenv )
TYPE ( global_environment_type ), INTENT ( IN ) :: globenv
INTEGER :: iw, ierror, ilen
CHARACTER ( LEN = 20 ) :: string
CHARACTER ( LEN = 5 ) :: label
!
iw = globenv % scr
!..defaults
runtest = 0
max_memory = 256.e6_dbl
!..parse the input section
label = '&TEST'
CALL parser_init ( globenv % input_file_name, label, ierror, globenv )
IF (ierror /= 0 ) THEN
IF( globenv % ionode ) &
WRITE ( iw, '( a )' ) ' No input section &TEST found '
ELSE
CALL read_line
DO WHILE ( test_next() /= 'X' )
ilen = 6
CALL cfield ( string, ilen )
CALL uppercase ( string )
SELECT CASE ( string )
CASE DEFAULT
CALL p_error()
CALL stop_parser ( 'test_input', 'unknown option' )
CASE ( 'MEMORY' )
max_memory = get_real ()
CASE ( 'COPY' )
runtest ( 1 ) = get_int ()
CASE ( 'MATMUL' )
runtest ( 2 ) = get_int ()
CASE ( 'FFT' )
runtest ( 3 ) = get_int ()
END SELECT
CALL read_line
END DO
END IF
CALL parser_end
!..end of parsing the input section
END SUBROUTINE test_input
!!*****
!******************************************************************************
!!****** library_tests/copy_test [1.0] *
!!
!! NAME
!! copy_test
!!
!! FUNCTION
!! Tests the performance to copy two vectors.
!!
!! AUTHOR
!! JGH 6-NOV-2000
!!
!! MODIFICATION HISTORY
!! none
!!
!! NOTES
!! The results of these tests allow to determine the size of the cache
!! of the CPU. This can be used to optimize the performance of the
!! FFTSG library.
!!
!! SOURCE
!******************************************************************************
SUBROUTINE copy_test ( globenv )
@ -57,14 +191,16 @@ SUBROUTINE copy_test ( globenv )
TYPE ( global_environment_type ), INTENT ( IN ) :: globenv
REAL ( dbl ), DIMENSION ( : ), ALLOCATABLE :: ca, cb
INTEGER :: len, ntim, iw, ierr, i, j
INTEGER :: len, ntim, iw, ierr, i, j, siz
REAL ( dbl ) :: perf, tstart, tend, t
! test for copy --> Cache size
siz = ABS ( runtest ( 1 ) )
iw = globenv % scr
IF ( globenv % ionode ) WRITE ( iw, '(//,A,/)' ) " Test of copy ( F95 ) "
DO i = 6, 24
len = 2**i
IF ( 8._dbl * REAL ( len, dbl ) > max_memory * 0.5_dbl ) EXIT
ALLOCATE ( ca ( len ), STAT = ierr )
IF ( ierr /= 0 ) EXIT
ALLOCATE ( cb ( len ), STAT = ierr )
@ -73,7 +209,7 @@ SUBROUTINE copy_test ( globenv )
CALL random_number ( ca )
ntim = NINT ( 1.e7_dbl / REAL ( len, dbl ) )
ntim = MAX ( ntim, 1 )
ntim = MIN ( ntim, 50000 )
ntim = MIN ( ntim, siz * 10000 )
tstart = m_cputime ( )
DO j = 1, ntim
@ -97,6 +233,24 @@ SUBROUTINE copy_test ( globenv )
END SUBROUTINE copy_test
!!*****
!******************************************************************************
!!****** library_tests/matmul_test [1.0] *
!!
!! NAME
!! matmul_test
!!
!! FUNCTION
!! Tests the performance of different kinds of matrix matrix multiply
!! kernels for the BLAS and F95 intrinsic matmul.
!!
!! AUTHOR
!! JGH 6-NOV-2000
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
SUBROUTINE matmul_test ( globenv )
@ -104,14 +258,16 @@ SUBROUTINE matmul_test ( globenv )
TYPE ( global_environment_type ), INTENT ( IN ) :: globenv
REAL ( dbl ), DIMENSION ( :, : ), ALLOCATABLE :: ma, mb, mc
INTEGER :: len, ntim, iw, ierr, i, j
INTEGER :: len, ntim, iw, ierr, i, j, siz
REAL ( dbl ) :: perf, tstart, tend, t
! test for matrix multpies
siz = ABS ( runtest ( 2 ) )
iw = globenv % scr
IF ( globenv % ionode ) WRITE ( iw, '(//,A,/)' ) " Test of matmul ( F95 ) "
DO i = 4, 10, 2
len = 2**i + 1
IF ( 8._dbl * REAL ( len*len, dbl ) > max_memory * 0.3_dbl ) EXIT
ALLOCATE ( ma ( len, len ), STAT = ierr )
IF ( ierr /= 0 ) EXIT
ALLOCATE ( mb ( len, len ), STAT = ierr )
@ -123,7 +279,7 @@ SUBROUTINE matmul_test ( globenv )
CALL random_number ( mb )
ntim = NINT ( 1.e8_dbl / ( 2._dbl * REAL ( len, dbl )**3 ) )
ntim = MAX ( ntim, 1 )
ntim = MIN ( ntim, 1000 )
ntim = MIN ( ntim, siz * 200 )
tstart = m_cputime ( )
DO j = 1, ntim
@ -193,6 +349,7 @@ SUBROUTINE matmul_test ( globenv )
IF ( globenv % ionode ) WRITE ( iw, '(//,A,/)' ) " Test of matmul ( BLAS ) "
DO i = 4, 10, 2
len = 2**i + 1
IF ( 8._dbl * REAL ( len*len, dbl ) > max_memory * 0.3_dbl ) EXIT
ALLOCATE ( ma ( len, len ), STAT = ierr )
IF ( ierr /= 0 ) EXIT
ALLOCATE ( mb ( len, len ), STAT = ierr )
@ -269,6 +426,23 @@ SUBROUTINE matmul_test ( globenv )
END SUBROUTINE matmul_test
!!*****
!******************************************************************************
!!****** library_tests/fft_test [1.0] *
!!
!! NAME
!! fft_test
!!
!! FUNCTION
!! Tests the performance of all available FFT libraries for 3D FFTs
!!
!! AUTHOR
!! JGH 6-NOV-2000
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
SUBROUTINE fft_test ( globenv )
@ -280,13 +454,14 @@ SUBROUTINE fft_test ( globenv )
COMPLEX ( dbl ), DIMENSION ( 4, 4, 4 ) :: zz
INTEGER :: len, ntim, iw, ierr, i, j, it, n(3), iall
REAL ( dbl ) :: perf, tstart, tend, t, flops, scale
INTEGER :: radix_in, radix_out, stat
INTEGER :: ndate ( 3 ) = (/ 12, 48, 128 /)
CHARACTER ( LEN = 6 ) :: method
INTEGER :: radix_in, radix_out, stat, siz
INTEGER :: ndate ( 3 ) = (/ 12, 48, 96 /)
CHARACTER ( LEN = 7 ) :: method
! test for 3d FFT
iw = globenv % scr
IF ( globenv % ionode ) WRITE ( iw, '(//,A,/)' ) " Test of 3D-FFT "
siz = ABS ( runtest ( 3 ) )
DO iall = 1, 100
SELECT CASE ( iall )
@ -294,10 +469,13 @@ SUBROUTINE fft_test ( globenv )
EXIT
CASE ( 1 )
CALL init_fft ( "FFTSG" )
method = "FFTSG "
method = "FFTSG "
CASE ( 2 )
CALL init_fft ( "FFTW" )
method = "FFTW "
method = "FFTW "
CASE ( 3 )
CALL init_fft ( "FFTESSL" )
method = "FFTESSL"
END SELECT
n = 4
CALL fft3d ( 1, n, zz, status=stat )
@ -307,6 +485,7 @@ SUBROUTINE fft_test ( globenv )
CALL fft_radix_operations ( radix_in, radix_out, FFT_RADIX_CLOSEST )
len = radix_out
n = len
IF ( 16._dbl * REAL ( len*len*len, dbl ) > max_memory * 0.5_dbl ) EXIT
ALLOCATE ( ra ( len, len, len ), STAT = ierr )
IF ( ierr /= 0 ) CALL stop_memory ( "lib_test", "ra", len*len*len )
ALLOCATE ( ca ( len, len, len ), STAT = ierr )
@ -316,9 +495,9 @@ SUBROUTINE fft_test ( globenv )
CALL random_number ( ra )
ca = ca + CMPLX ( 0._dbl, 1._dbl ) * ra
flops = REAL ( len**3, dbl ) * 15._dbl * LOG ( REAL ( len, dbl ) )
ntim = NINT ( 5.e7_dbl / flops )
ntim = NINT ( siz * 1.e7_dbl / flops )
ntim = MAX ( ntim, 1 )
ntim = MIN ( ntim, 1000 )
ntim = MIN ( ntim, 200 )
scale = 1._dbl / REAL ( len**3, dbl )
tstart = m_cputime ( )
DO j = 1, ntim
@ -345,6 +524,7 @@ SUBROUTINE fft_test ( globenv )
END SUBROUTINE fft_test
!!*****
!******************************************************************************
END MODULE library_tests

18
src/linklist_cell_list.F
View file

@ -2,6 +2,21 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/linklist_cell_list [1.0] *
!!
!! NAME
!! linklist_cell_list
!!
!! FUNCTION
!!
!! AUTHOR
!! CJM
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE linklist_cell_list
@ -20,6 +35,9 @@ MODULE linklist_cell_list
PRIVATE
PUBLIC :: cell_list
!!*****
!******************************************************************************
CONTAINS
!******************************************************************************

18
src/linklist_cell_types.F
View file

@ -2,6 +2,21 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****s* cp2k/linklist_cell_types [1.0] *
!!
!! NAME
!! linklist_cell_types
!!
!! FUNCTION
!!
!! AUTHOR
!! CJM
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE linklist_cell_types
@ -34,3 +49,6 @@ MODULE linklist_cell_types
END TYPE cell_neighborlist_type
END MODULE linklist_cell_types
!******************************************************************************
!!*****

162
src/linklist_control.F
View file

@ -2,6 +2,22 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2k developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/linklist_control [1.0] *
!!
!! NAME
!! linklist_control
!!
!! FUNCTION
!!
!! AUTHOR
!! CJM
!!
!! MODIFICATION HISTORY
!! Harald Forbert (Dec-2000): Changes for multiple linked lists
!! linklist_internal_data_type
!!
!! SOURCE
!******************************************************************************
MODULE linklist_control
@ -20,47 +36,28 @@ MODULE linklist_control
USE simulation_cell, ONLY : cell_type
USE stop_program, ONLY : stop_prg, stop_memory
USE timings, ONLY : timeset, timestop
USE linklist_types, ONLY : linklist_internal_data_type
USE structure_types, ONLY : structure_type
IMPLICIT NONE
PRIVATE
PUBLIC :: set_ll_parm, get_ll_parm, list_control
! globals specific to the particle-particle verlet list
TYPE ( cell_neighborlist_type ), DIMENSION (:,:,:), POINTER :: pp_cell_ll
TYPE ( cell_atoms_linklist_type ), DIMENSION (:,:,:,:), POINTER :: &
pp_startlist_cell
TYPE ( neighbor_cells_type ), DIMENSION ( :, : ), POINTER :: pp_neighbor
TYPE ( linklist_neighbor ), POINTER :: pp_startlist_nl
TYPE ( startlist_images_type ), POINTER :: pp_startlist_im ( : )
INTEGER, DIMENSION ( 3 ) :: pp_cell
INTEGER, POINTER :: pp_images(:,:,:)
INTEGER, POINTER :: pp_ncell ( : )
TYPE ( global_environment_type ) :: globenv_ll
INTEGER :: natom_types=0
INTEGER :: list_type = 0
INTEGER :: counter = 0, last_update = 0, num_update = 0
INTEGER :: print_level = 1
REAL ( dbl ), DIMENSION ( :, : ), ALLOCATABLE :: rlist_cut, rlist_cutsq
REAL ( dbl ) :: verlet_skin
REAL ( dbl ) :: lup, aup
!print_level > 0 : no output
!print_level < 5 : normal output
!print_level < 10 : initial neighbor list info
!print_level > 10 : full neighbor list info at each update
!!*****
!******************************************************************************
CONTAINS
!******************************************************************************
SUBROUTINE set_ll_parm ( globenv_in, vskin, natype, rcut, &
SUBROUTINE set_ll_parm ( internal_data, globenv_in, vskin, natype, rcut, &
ppcells, count, printlevel, ltype )
IMPLICIT NONE
! Arguments
TYPE ( linklist_internal_data_type ), INTENT ( INOUT ) :: internal_data
TYPE ( global_environment_type ), INTENT ( IN ) :: globenv_in
REAL ( dbl ), INTENT ( IN ), OPTIONAL :: vskin
REAL ( dbl ), INTENT ( IN ), OPTIONAL :: rcut ( :, : )
@ -77,52 +74,52 @@ SUBROUTINE set_ll_parm ( globenv_in, vskin, natype, rcut, &
!------------------------------------------------------------------------------
globenv_ll = globenv_in
internal_data%globenv_ll = globenv_in
! Defines the number of cells for the particle-particle lists
IF ( PRESENT ( ppcells)) pp_cell = ppcells
IF ( PRESENT ( ppcells)) internal_data%pp_cell = ppcells
! Defines the verlet skin
IF ( PRESENT ( vskin)) verlet_skin = vskin
IF ( PRESENT ( vskin)) internal_data%verlet_skin = vskin
! Defines the number of atom types
IF ( PRESENT ( natype )) natom_types = natype
IF ( PRESENT ( natype )) internal_data%natom_types = natype
! Defines cut-off for particle-particle lists
IF ( PRESENT ( rcut)) THEN
! allocating rlist and rlist_cut to upper triangular part
! of interaction matrix
ALLOCATE (rlist_cut(natom_types,natom_types),STAT=isos)
ALLOCATE (internal_data%rlist_cut(internal_data%natom_types,internal_data%natom_types),STAT=isos)
IF ( isos /= 0 ) CALL stop_memory ( 'set_ll_param', &
'rlist_cut', natom_types ** 2 )
'rlist_cut', internal_data%natom_types ** 2 )
ALLOCATE ( rlist_cutsq ( natom_types, natom_types ), STAT = isos )
ALLOCATE ( internal_data%rlist_cutsq ( internal_data%natom_types, internal_data%natom_types ), STAT = isos )
IF ( isos /= 0 ) CALL stop_memory ( 'set_ll_param', &
'rlist_cutsq', natom_types ** 2 )
DO iat = 1, natom_types
DO jat = iat, natom_types
rlist_cut(iat,jat) = sqrt(rcut(iat,jat)) + verlet_skin
rlist_cutsq(iat,jat) = rlist_cut(iat,jat)*rlist_cut(iat,jat)
rlist_cut(jat,iat) = rlist_cut(iat,jat)
rlist_cutsq(jat,iat) = rlist_cutsq(iat,jat)
'rlist_cutsq', internal_data%natom_types ** 2 )
DO iat = 1, internal_data%natom_types
DO jat = iat, internal_data%natom_types
internal_data%rlist_cut(iat,jat) = sqrt(rcut(iat,jat)) + internal_data%verlet_skin
internal_data%rlist_cutsq(iat,jat) = internal_data%rlist_cut(iat,jat)*internal_data%rlist_cut(iat,jat)
internal_data%rlist_cut(jat,iat) = internal_data%rlist_cut(iat,jat)
internal_data%rlist_cutsq(jat,iat) = internal_data%rlist_cutsq(iat,jat)
END DO
END DO
END IF
! Defines the counting variable
IF ( PRESENT ( count ) ) counter = count
IF ( PRESENT ( count ) ) internal_data%counter = count
! Defines the print level
IF ( PRESENT ( printlevel ) ) print_level = printlevel
IF ( PRESENT ( printlevel ) ) internal_data%print_level = printlevel
! Defines the type of list to be constructed
IF ( PRESENT ( ltype ) ) THEN
SELECT CASE ( ltype )
CASE DEFAULT
list_type = 0
internal_data%list_type = 0
CASE ( 'NONBOND')
list_type = 1
internal_data%list_type = 1
END SELECT
END IF
@ -130,14 +127,15 @@ END SUBROUTINE set_ll_parm
!******************************************************************************
SUBROUTINE get_ll_parm ( last_up, av_up )
SUBROUTINE get_ll_parm ( internal_data, last_up, av_up )
IMPLICIT NONE
TYPE ( linklist_internal_data_type ), INTENT ( IN ) :: internal_data
REAL ( dbl ), INTENT ( OUT ) :: last_up, av_up
last_up = lup
av_up = aup
last_up = internal_data%lup
av_up = internal_data%aup
END SUBROUTINE get_ll_parm
@ -149,18 +147,18 @@ 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 ( pnode, part, box )
SUBROUTINE list_control ( internal_data, pnode, part, box )
IMPLICIT NONE
! Arguments
TYPE ( linklist_internal_data_type ), INTENT ( INOUT ) :: internal_data
TYPE ( particle_type ), DIMENSION ( : ), INTENT ( INOUT ) :: part
TYPE ( particle_node_type ), DIMENSION ( : ), INTENT ( INOUT ) :: pnode
TYPE ( cell_type ), INTENT ( IN ) :: box
! Locals
REAL ( dbl ), DIMENSION (3) :: disp_vec
REAL ( dbl ), ALLOCATABLE, SAVE :: r_last_update ( :, : )
REAL ( dbl ) :: displace, max_displace
LOGICAL :: list_update_flag, ionode
INTEGER :: i, nnodes, isos, iw, handle
@ -171,38 +169,38 @@ SUBROUTINE list_control ( pnode, part, box )
CALL timeset ( 'LIST_CONTROL','I',' ',handle )
nnodes = SIZE ( pnode )
IF ( .NOT. ALLOCATED ( r_last_update ) ) THEN
ALLOCATE ( r_last_update ( 3, nnodes ), STAT = isos )
IF ( .NOT. ASSOCIATED ( internal_data%r_last_update ) ) THEN
ALLOCATE ( internal_data%r_last_update ( 3, nnodes ), STAT = isos )
IF ( isos /= 0 ) CALL stop_memory ( 'list_control', &
'r_last_update', 3 * nnodes )
r_last_update = 1.0E+10_dbl
internal_data%r_last_update = 1.0E+10_dbl
END IF
! check to see if neighbor list has to be updated
max_displace = -1.0E+15_dbl
DO i = 1, nnodes
disp_vec ( : ) = r_last_update(:,i) - pnode(i) %p%r ( : )
disp_vec ( : ) = internal_data%r_last_update(:,i) - pnode(i) %p%r ( : )
displace = DOT_PRODUCT ( disp_vec, disp_vec )
max_displace = MAX ( displace, max_displace )
END DO
#if defined(__parallel)
CALL mp_max ( max_displace, globenv_ll%group )
CALL mp_max ( max_displace, internal_data%globenv_ll%group )
#endif
IF ( max_displace >= 0.25_dbl * verlet_skin ** 2 ) THEN
IF ( max_displace >= 0.25_dbl * internal_data%verlet_skin ** 2 ) THEN
list_update_flag = .TRUE.
DO i = 1, nnodes
r_last_update(1,i) = pnode(i) % p % r(1)
r_last_update(2,i) = pnode(i) % p % r(2)
r_last_update(3,i) = pnode(i) % p % r(3)
internal_data%r_last_update(1,i) = pnode(i) % p % r(1)
internal_data%r_last_update(2,i) = pnode(i) % p % r(2)
internal_data%r_last_update(3,i) = pnode(i) % p % r(3)
END DO
ELSE
list_update_flag = .FALSE.
END IF
IF ( list_update_flag ) THEN
SELECT CASE ( list_type )
SELECT CASE ( internal_data%list_type )
CASE ( 0 )
CALL stop_prg ( 'list_control', 'list type not specified' )
@ -210,33 +208,39 @@ SUBROUTINE list_control ( pnode, part, box )
CASE ( 1 )
string = ' cell_index: nonbonded neighbor list'
CALL cell_list(natom_types,pp_cell,rlist_cut,string,pnode,part, &
box,pp_cell_ll,pp_startlist_cell,pp_neighbor, globenv_ll)
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)
CALL verlet_list(natom_types,pp_images,pp_cell,pnode,part,box, &
verlet_skin,rlist_cut,rlist_cutsq,pp_neighbor,pp_cell_ll, &
pp_startlist_nl,pp_startlist_im, globenv_ll)
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)
END SELECT
iw = globenv_ll % scr
ionode = globenv_ll % ionode
IF ( ionode .AND. print_level > 5 ) WRITE (iw,'(a,T71,I10 )' ) &
' list_control| Neighbor list updated at step: ', counter
IF ( counter > 0 ) THEN
num_update = num_update + 1
lup = counter + 1 - last_update
last_update = counter + 1
aup = aup + (lup-aup)/float(num_update)
iw = internal_data%globenv_ll % scr
ionode = internal_data%globenv_ll % ionode
IF ( ionode .AND. internal_data%print_level > 5 ) WRITE (iw,'(a,T71,I10 )' ) &
' list_control| Neighbor list updated at step: ', internal_data%counter
IF ( internal_data%counter > 0 ) THEN
internal_data%num_update = internal_data%num_update + 1
internal_data%lup = internal_data%counter + 1 - internal_data%last_update
internal_data%last_update = internal_data%counter + 1
internal_data%aup = internal_data%aup + &
(internal_data%lup-internal_data%aup)/float(internal_data%num_update)
ELSE
num_update = 0
lup = 0
last_update = 1
aup = 0.0_dbl
internal_data%num_update = 0
internal_data%lup = 0
internal_data%last_update = 1
internal_data%aup = 0.0_dbl
END IF
END IF
counter = counter + 1
internal_data%counter = internal_data%counter + 1
CALL timestop ( 0._dbl, handle )
@ -245,3 +249,5 @@ END SUBROUTINE list_control
!******************************************************************************
END MODULE linklist_control
!******************************************************************************

61
src/linklist_types.F Normal file
View file

@ -0,0 +1,61 @@
!-----------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****s* cp2k/linklist_types [1.0] *
!!
!! NAME
!! linklist_types
!!
!! FUNCTION
!!
!! AUTHOR
!! HAF
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE linklist_types
USE linklist_cell_types, ONLY : cell_neighborlist_type, &
cell_atoms_linklist_type, &
neighbor_cells_type
USE molecule_types, ONLY : linklist_neighbor
USE linklist_utilities, ONLY : startlist_images_type
USE global_types, ONLY : global_environment_type
USE kinds, ONLY : dbl
IMPLICIT NONE
PRIVATE
PUBLIC :: linklist_internal_data_type
TYPE linklist_internal_data_type
TYPE ( cell_neighborlist_type ), DIMENSION (:,:,:), POINTER :: pp_cell_ll
TYPE ( cell_atoms_linklist_type ), DIMENSION (:,:,:,:), POINTER :: &
pp_startlist_cell
TYPE ( neighbor_cells_type ), DIMENSION ( :, : ), POINTER :: pp_neighbor
TYPE ( linklist_neighbor ), POINTER :: pp_startlist_nl
TYPE ( startlist_images_type ), POINTER :: pp_startlist_im ( : )
INTEGER, DIMENSION ( 3 ) :: pp_cell
INTEGER, POINTER :: pp_images(:,:,:)
INTEGER, POINTER :: pp_ncell ( : )
TYPE ( global_environment_type ) :: globenv_ll
INTEGER :: natom_types
INTEGER :: list_type
INTEGER :: counter, last_update, num_update
INTEGER :: print_level
REAL ( dbl ), DIMENSION ( :, : ), POINTER :: rlist_cut
REAL ( dbl ), DIMENSION ( :, : ), POINTER :: rlist_cutsq
REAL ( dbl ) :: verlet_skin
REAL ( dbl ) :: lup, aup
REAL ( dbl ), DIMENSION ( :, : ), POINTER :: r_last_update
END TYPE linklist_internal_data_type
END MODULE linklist_types
!******************************************************************************
!!*****

17
src/linklist_utilities.F
View file

@ -2,6 +2,21 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/linklist_utilities [1.0] *
!!
!! NAME
!! linklist_utilities
!!
!! FUNCTION
!!
!! AUTHOR
!! CJM
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE linklist_utilities
@ -22,6 +37,8 @@ MODULE linklist_utilities
TYPE ( linklist_images), POINTER :: image
END TYPE startlist_images_type
!!*****
!------------------------------------------------------------------------------!
CONTAINS

18
src/linklist_verlet_list.F
View file

@ -2,6 +2,21 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/linklist_verlet_list [1.0] *
!!
!! NAME
!! linklist_verlet_list
!!
!! FUNCTION
!!
!! AUTHOR
!! CJM
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE linklist_verlet_list
@ -25,6 +40,9 @@ MODULE linklist_verlet_list
PRIVATE
PUBLIC :: verlet_list
!!*****
!******************************************************************************
CONTAINS
!******************************************************************************

142
src/linklists.F
View file

@ -2,6 +2,22 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/linklists [1.0] *
!!
!! NAME
!! linklists
!!
!! FUNCTION
!! Generate linked lists for intramolecular dependencies
!!
!! AUTHOR
!! CJM & JGH
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE linklists
@ -22,9 +38,28 @@ MODULE linklists
PUBLIC :: bonds, bends, torsions, dist_constraints, g3x3_constraints
LOGICAL :: bonds, bends, torsions, dist_constraints, g3x3_constraints
!!*****
!******************************************************************************
CONTAINS
!******************************************************************************
!!****** linklists/atom_link_list [1.0] *
!!
!! NAME
!! atom_link_list
!!
!! FUNCTION
!! Generate linked lists for intramolecular dependencies
!!
!! AUTHOR
!! CJM & JGH
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
SUBROUTINE atom_link_list ( molecule, pnode )
@ -61,6 +96,23 @@ SUBROUTINE atom_link_list ( molecule, pnode )
END SUBROUTINE atom_link_list
!!*****
!******************************************************************************
!!****** linklists/exclusion [1.0] *
!!
!! NAME
!! exclusion
!!
!! FUNCTION
!! Generate linked lists for intramolecular dependencies
!!
!! AUTHOR
!! CJM & JGH
!!
!! MODIFICATION HISTORY
!! Debugged parallel code, (JGH) 25-Nov-2000
!!
!! SOURCE
!******************************************************************************
SUBROUTINE exclusion(molecule,pnode)
@ -174,7 +226,7 @@ SUBROUTINE exclusion(molecule,pnode)
END DO
BENDIF: IF (match) THEN
DO ii = 1, 3
CALL checklist(pnode(jatom),match,jatom,bend(ii))
CALL checklist(pnode(jnode),match,jatom,bend(ii))
IF ( .NOT. match) THEN
i = getnode(molecule,bend(ii))
current_excl%p => pnode(i) %p
@ -366,6 +418,23 @@ FUNCTION getnode ( molecule, atom ) RESULT ( node )
END FUNCTION getnode
!!*****
!******************************************************************************
!!****** linklists/bond_link_list [1.0] *
!!
!! NAME
!! bond_link_list
!!
!! FUNCTION
!! Generate linked lists for intramolecular dependencies
!!
!! AUTHOR
!! CJM & JGH
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
SUBROUTINE bond_link_list ( pnode, molecule, mol_setup, bond_parm )
@ -425,6 +494,23 @@ SUBROUTINE bond_link_list ( pnode, molecule, mol_setup, bond_parm )
END SUBROUTINE bond_link_list
!!*****
!******************************************************************************
!!****** linklists/bend_link_list [1.0] *
!!
!! NAME
!! bend_link_list
!!
!! FUNCTION
!! Generate linked lists for intramolecular dependencies
!!
!! AUTHOR
!! CJM & JGH
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
SUBROUTINE bend_link_list ( pnode, molecule, mol_setup, bend_parm )
@ -488,6 +574,23 @@ SUBROUTINE bend_link_list ( pnode, molecule, mol_setup, bend_parm )
END SUBROUTINE bend_link_list
!!*****
!******************************************************************************
!!****** linklists/torsion_link_list [1.0] *
!!
!! NAME
!! torsion_link_list
!!
!! FUNCTION
!! Generate linked lists for intramolecular dependencies
!!
!! AUTHOR
!! CJM & JGH
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
SUBROUTINE torsion_link_list ( pnode, molecule, mol_setup, torsion_parm )
@ -557,6 +660,23 @@ SUBROUTINE torsion_link_list ( pnode, molecule, mol_setup, torsion_parm )
END SUBROUTINE torsion_link_list
!!*****
!******************************************************************************
!!****** linklists/distconst_link_list [1.0] *
!!
!! NAME
!! distconst_link_list
!!
!! FUNCTION
!! Generate linked lists for intramolecular dependencies
!!
!! AUTHOR
!! CJM & JGH
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
SUBROUTINE distconst_link_list ( pnode, molecule, mol_setup, &
@ -621,6 +741,23 @@ SUBROUTINE distconst_link_list ( pnode, molecule, mol_setup, &
END SUBROUTINE distconst_link_list
!!*****
!******************************************************************************
!!****** linklists/g3x3_link_list [1.0] *
!!
!! NAME
!! g3x3_link_list
!!
!! FUNCTION
!! Generate linked lists for intramolecular dependencies
!!
!! AUTHOR
!! CJM & JGH
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
SUBROUTINE g3x3_link_list ( pnode, molecule, mol_setup, constraint_distance )
@ -688,6 +825,9 @@ SUBROUTINE g3x3_link_list ( pnode, molecule, mol_setup, constraint_distance )
END SUBROUTINE g3x3_link_list
!!*****
!******************************************************************************
END MODULE linklists
!******************************************************************************

17
src/machine.F
View file

@ -2,6 +2,21 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/machine [1.0] *
!!
!! NAME
!! machine
!!
!! FUNCTION
!!
!! AUTHOR
!! APSI & JGH
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE machine
@ -44,3 +59,5 @@ MODULE machine
!
END MODULE machine
!!*****

20
src/machine_absoft.F
View file

@ -2,6 +2,23 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/machine_absoft [1.0] *
!!
!! NAME
!! machine_absoft
!!
!! FUNCTION
!!
!! AUTHOR
!! APSI & JGH
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!!*****
!******************************************************************************
MODULE machine_absoft
USE kinds, ONLY : dbl, sgl
@ -10,6 +27,9 @@ MODULE machine_absoft
PUBLIC :: m_walltime, m_cputime, m_datum
PUBLIC :: m_hostnm, m_getcwd, m_getlog, m_getuid, m_getpid, m_getarg
!!*****
!******************************************************************************
CONTAINS
!******************************************************************************

18
src/machine_aix.F
View file

@ -2,6 +2,21 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/machine_aix [1.0] *
!!
!! NAME
!! machine_aix
!!
!! FUNCTION
!!
!! AUTHOR
!! APSI & JGH
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE machine_aix
@ -11,6 +26,9 @@ MODULE machine_aix
PUBLIC :: m_walltime, m_cputime, m_datum, print_memory
PUBLIC :: m_hostnm, m_getcwd, m_getlog, m_getuid, m_getpid, m_getarg
!!*****
!******************************************************************************
CONTAINS
!******************************************************************************

18
src/machine_dec.F
View file

@ -2,6 +2,21 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/machine_dec [1.0] *
!!
!! NAME
!! machine_dec
!!
!! FUNCTION
!!
!! AUTHOR
!! APSI & JGH
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE machine_dec
@ -11,6 +26,9 @@ MODULE machine_dec
PUBLIC :: m_walltime, m_cputime, m_datum, print_memory
PUBLIC :: m_hostnm, m_getcwd, m_getlog, m_getuid, m_getpid, m_getarg
!!*****
!******************************************************************************
CONTAINS
!******************************************************************************

18
src/machine_irix.F
View file

@ -2,6 +2,21 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/machine_irix [1.0] *
!!
!! NAME
!! machine_irix
!!
!! FUNCTION
!!
!! AUTHOR
!! APSI & JGH
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE machine_irix
@ -11,6 +26,9 @@ MODULE machine_irix
PUBLIC :: m_walltime, m_cputime, m_datum
PUBLIC :: m_hostnm, m_getcwd, m_getlog, m_getuid, m_getpid, m_getarg
!!*****
!******************************************************************************
CONTAINS
!******************************************************************************

19
src/machine_pgi.F
View file

@ -2,6 +2,22 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/machine_pgi [1.0] *
!!
!! NAME
!! machine_pgi
!!
!! FUNCTION
!!
!! AUTHOR
!! APSI & JGH
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!!*****
!******************************************************************************
MODULE machine_pgi
@ -11,6 +27,9 @@ MODULE machine_pgi
PUBLIC :: m_walltime, m_cputime, m_datum, print_memory
PUBLIC :: m_hostnm, m_getcwd, m_getlog, m_getuid, m_getpid, m_getarg
!!*****
!******************************************************************************
CONTAINS
!******************************************************************************

18
src/machine_t3e.F
View file

@ -2,6 +2,21 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/machine_t3e [1.0] *
!!
!! NAME
!! machine_t3e
!!
!! FUNCTION
!!
!! AUTHOR
!! APSI & JGH
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE machine_t3e
@ -11,6 +26,9 @@ MODULE machine_t3e
PUBLIC :: m_walltime, m_cputime, m_datum, print_memory
PUBLIC :: m_hostnm, m_getcwd, m_getlog, m_getuid, m_getpid, m_getarg
!!*****
!******************************************************************************
CONTAINS
!******************************************************************************

39
src/mathconstants.F
View file

@ -2,7 +2,22 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!! Definition of mathematical constants and functions.
!!****** cp2k/mathconstants [1.0] *
!!
!! NAME
!! mathconstants
!!
!! FUNCTION
!! Definition of mathematical constants and functions.
!!
!! AUTHOR
!! Matthias Krack
!!
!! MODIFICATION HISTORY
!! Adapted for use in CP2K (JGH)
!!
!! SOURCE
!******************************************************************************
MODULE mathconstants
@ -32,8 +47,27 @@ MODULE mathconstants
REAL ( dbl ), DIMENSION (0:maxfac) :: fac
REAL ( dbl ), DIMENSION (-1:2*maxfac+1) :: dfac
!!*****
!******************************************************************************
CONTAINS
!******************************************************************************
!!****** cp2k/init_mathcon [1.0] *
!!
!! NAME
!! init_mathcon
!!
!! FUNCTION
!! Definition of mathematical constants and functions.
!!
!! AUTHOR
!! JGH
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
SUBROUTINE init_mathcon
@ -66,6 +100,9 @@ SUBROUTINE init_mathcon
END SUBROUTINE init_mathcon
!!*****
!******************************************************************************
END MODULE mathconstants
!******************************************************************************

137
src/md.F
View file

@ -2,6 +2,21 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/md [1.0] *
!!
!! NAME
!! md
!!
!! FUNCTION
!!
!! AUTHOR
!! CJM
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE md
@ -56,48 +71,68 @@ MODULE md
CHARACTER ( LEN = 50 ) :: dump_file_name
END TYPE mdio_parameters_type
!!*****
!******************************************************************************
CONTAINS
!!>---------------------------------------------------------------------------!
!! SECTION: &md ... &end !
!! !
!! ensemble: [nve,nvt] ensemble type !
!! steps: n total number of steps !
!! timestep: dt time step [fs] !
!! temperature: T temperature [K] !
!! restart: [pos,posvel] start type !
!! constraints: [on,off] constraints !
!! shake: eps convergence for Shake algorithm !
!! verlet_skin vs buffer size for Verlet lists !
!! cells x y z Number of cells in each direction !
!! Nose_parameter & parameter of Nose-Hoover-chains !
!! length: nhclen & length of thermostat !
!! Yoshida: nyosh & order of Yoshida integrator !
!! mts: nc & number of multiple time steps !
!! timecon: tau_nhc NHC time constant !
!! Barostat_parameter parameter of barostat !
!! pressure: P & pressure [bar] !
!! timecon: tau_cell & barostat time constant !
!! FILES & !
!! filebody "filename" & !
!! coordinates "filename.crd" & !
!! velocities "filename.vel" & !
!! pressure "filename.ptens" & !
!! energies "filename.ener" & !
!! temperature "filename.temp" & !
!! dump "filename.dum" !
!! PRINT & !
!! screen iscreen & !
!! files ifiles & !
!! coordinates icrd & !
!! velocities ivel & !
!! pressure iptens & !
!! energies iener & !
!! temperature itemp & !
!! dump idump !
!! !
!! !
!!<---------------------------------------------------------------------------!
!******************************************************************************
!!****** md/read_md_section [1.0] *
!!
!! NAME
!! read_md_section
!!
!! FUNCTION
!!
!! AUTHOR
!! CJM
!!
!! MODIFICATION HISTORY
!! none
!!
!! NOTES
!!I---------------------------------------------------------------------------I
!!I SECTION: &md ... &end I
!!I I
!!I 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 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 Nose_parameter & parameter of Nose-Hoover-chains 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 Barostat_parameter parameter of barostat I
!!I pressure: P & pressure [bar] I
!!I timecon: tau_cell & barostat time constant I
!!I FILES & I
!!I filebody "filename" & I
!!I coordinates "filename.crd" & I
!!I velocities "filename.vel" & I
!!I pressure "filename.ptens" & I
!!I energies "filename.ener" & I
!!I temperature "filename.temp" & I
!!I dump "filename.dum" I
!!I PRINT & I
!!I screen iscreen & I
!!I files ifiles & I
!!I coordinates icrd & I
!!I velocities ivel & I
!!I pressure iptens & I
!!I energies iener & I
!!I temperature itemp & I
!!I dump idump I
!!I I
!!I I
!!I---------------------------------------------------------------------------I
!! SOURCE
!******************************************************************************
SUBROUTINE read_md_section ( simpar, mdpar, mdio )
@ -397,9 +432,24 @@ SUBROUTINE read_md_section ( simpar, mdpar, mdio )
END SUBROUTINE read_md_section
!!*****
!******************************************************************************
!!****** md/initialize_velocities [1.0] *
!!
!! NAME
!! initialize_velocities
!!
!! FUNCTION
!! Initializes the velocities to the Maxwellian distribution
!!
!! AUTHOR
!! CJM
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
! Initializes the velocities to the Maxwellian distribution
SUBROUTINE initialize_velocities ( simpar, part, mdpar )
@ -494,6 +544,9 @@ SUBROUTINE initialize_velocities ( simpar, part, mdpar )
END SUBROUTINE initialize_velocities
!!*****
!******************************************************************************
END MODULE md
!******************************************************************************

18
src/message_passing.F
View file

@ -2,6 +2,21 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/message_passing [1.0] *
!!
!! NAME
!! message_passing
!!
!! FUNCTION
!!
!! AUTHOR
!! JGH
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE message_passing
@ -31,6 +46,9 @@ MODULE message_passing
MODULE PROCEDURE mp_min_i, mp_min_r, mp_min_rv
END INTERFACE
!!*****
!******************************************************************************
CONTAINS
!******************************************************************************

206
src/mol_force.F
View file

@ -2,6 +2,21 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/mol_force [1.0] *
!!
!! NAME
!! mol_force
!!
!! FUNCTION
!!
!! AUTHOR
!! CJM
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE mol_force
@ -11,6 +26,9 @@ MODULE mol_force
PUBLIC :: force_bonds, force_bends, force_torsions
PUBLIC :: get_pv_bond, get_pv_bend, get_pv_torsion
!!*****
!******************************************************************************
CONTAINS
!******************************************************************************
@ -74,108 +92,190 @@ END SUBROUTINE force_bends
!******************************************************************************
SUBROUTINE force_torsions
SUBROUTINE force_torsions ( is32, dist1, dist2, u1, u2, iw1, iw2, &
a, gt1, gt2, gt3, gt4, energy, fscalar )
IMPLICIT NONE
STOP 1991
! Arguments
REAL ( dbl ), INTENT ( IN ), DIMENSION ( : ) :: u1, u2
REAL ( dbl ), INTENT ( INOUT ), DIMENSION ( : ) :: gt1, gt2, gt3, gt4
REAL ( dbl ), INTENT ( IN ) :: is32, dist1, dist2, iw1, iw2
REAL ( dbl ), INTENT ( IN ), DIMENSION ( : ) :: a
REAL ( dbl ), INTENT ( OUT ) :: fscalar, energy
! Locals
REAL ( dbl ) :: cosphi
!------------------------------------------------------------------------------
! define
cosphi = DOT_PRODUCT ( u1, u2 ) * iw1 * iw2
! compute energy
energy = a(0) + a(1) + 2._dbl * a(3) &
+ ( a(1) - 3._dbl * a(3) ) * cosphi &
- 2._dbl * a(2) * cosphi ** 2 + 4._dbl * a(3) * cosphi ** 3
! compute fscalar
fscalar = - a(1) + 3._dbl * a(3) + 4._dbl * a(2) * cosphi - 12._dbl * a(3) * cosphi ** 2
! compute the gradients
gt1 = ( u2 * iw2 - cosphi * u1 * iw1 ) * iw1
gt4 = ( u1 * iw1 - cosphi * u2 * iw2 ) * iw2
gt2 = ( dist1 * is32 ** 2 - 1._dbl ) * gt1 + dist2 * is32 ** 2 * gt4
gt3 = ( - dist2 * is32 ** 2 - 1._dbl ) * gt4 - dist1 * is32 ** 2 * gt1
END SUBROUTINE force_torsions
!******************************************************************************
SUBROUTINE get_pv_bond(f1,f2,r1,r2,pvbond)
SUBROUTINE get_pv_bond(f1,f2,r1,r2,pv_bond)
IMPLICIT NONE
REAL ( dbl ), DIMENSION (3), INTENT ( IN ) :: f1, f2, r1, r2
REAL ( dbl ), DIMENSION ( :, : ), INTENT ( INOUT ) :: pvbond
REAL ( dbl ), DIMENSION ( :, : ), INTENT ( INOUT ) :: pv_bond
!------------------------------------------------------------------------------
! p(1,1)
pvbond(1,1) = pvbond(1,1) + f1(1)*r1(1)
pvbond(1,1) = pvbond(1,1) + f2(1)*r2(1)
pv_bond(1,1) = pv_bond(1,1) + f1(1)*r1(1)
pv_bond(1,1) = pv_bond(1,1) + f2(1)*r2(1)
! p(1,2)
pvbond(1,2) = pvbond(1,2) + f1(1)*r1(2)
pvbond(1,2) = pvbond(1,2) + f2(1)*r2(2)
pv_bond(1,2) = pv_bond(1,2) + f1(1)*r1(2)
pv_bond(1,2) = pv_bond(1,2) + f2(1)*r2(2)
! p(1,3)
pvbond(1,3) = pvbond(1,3) + f1(1)*r1(3)
pvbond(1,3) = pvbond(1,3) + f2(1)*r2(3)
pv_bond(1,3) = pv_bond(1,3) + f1(1)*r1(3)
pv_bond(1,3) = pv_bond(1,3) + f2(1)*r2(3)
! p(2,1)
pvbond(2,1) = pvbond(2,1) + f1(2)*r1(1)
pvbond(2,1) = pvbond(2,1) + f2(2)*r2(1)
pv_bond(2,1) = pv_bond(2,1) + f1(2)*r1(1)
pv_bond(2,1) = pv_bond(2,1) + f2(2)*r2(1)
! p(2,2)
pvbond(2,2) = pvbond(2,2) + f1(2)*r1(2)
pvbond(2,2) = pvbond(2,2) + f2(2)*r2(2)
pv_bond(2,2) = pv_bond(2,2) + f1(2)*r1(2)
pv_bond(2,2) = pv_bond(2,2) + f2(2)*r2(2)
! p(2,3)
pvbond(2,3) = pvbond(2,3) + f1(2)*r1(3)
pvbond(2,3) = pvbond(2,3) + f2(2)*r2(3)
pv_bond(2,3) = pv_bond(2,3) + f1(2)*r1(3)
pv_bond(2,3) = pv_bond(2,3) + f2(2)*r2(3)
! p(3,1)
pvbond(3,1) = pvbond(3,1) + f1(3)*r1(1)
pvbond(3,1) = pvbond(3,1) + f2(3)*r2(1)
pv_bond(3,1) = pv_bond(3,1) + f1(3)*r1(1)
pv_bond(3,1) = pv_bond(3,1) + f2(3)*r2(1)
! p(3,2)
pvbond(3,2) = pvbond(3,2) + f1(3)*r1(2)
pvbond(3,2) = pvbond(3,2) + f2(3)*r2(2)
pv_bond(3,2) = pv_bond(3,2) + f1(3)*r1(2)
pv_bond(3,2) = pv_bond(3,2) + f2(3)*r2(2)
! p(3,3)
pvbond(3,3) = pvbond(3,3) + f1(3)*r1(3)
pvbond(3,3) = pvbond(3,3) + f2(3)*r2(3)
pv_bond(3,3) = pv_bond(3,3) + f1(3)*r1(3)
pv_bond(3,3) = pv_bond(3,3) + f2(3)*r2(3)
END SUBROUTINE get_pv_bond
!******************************************************************************
SUBROUTINE get_pv_bend(f1,f2,f3,r1,r2,r3,pvbend)
SUBROUTINE get_pv_bend(f1,f2,f3,r1,r2,r3,pv_bend)
IMPLICIT NONE
! Arguments
REAL ( dbl ), DIMENSION (3), INTENT ( IN ) :: f1, f2, f3, r1, r2, r3
REAL ( dbl ), DIMENSION ( :, : ), INTENT ( INOUT ) :: pvbend
REAL ( dbl ), DIMENSION ( :, : ), INTENT ( INOUT ) :: pv_bend
!------------------------------------------------------------------------------
!p(1,1)
pvbend(1,1) = pvbend(1,1) + f1(1)*r1(1)
pvbend(1,1) = pvbend(1,1) + f2(1)*r2(1)
pvbend(1,1) = pvbend(1,1) + f3(1)*r3(1)
pv_bend(1,1) = pv_bend(1,1) + f1(1)*r1(1)
pv_bend(1,1) = pv_bend(1,1) + f2(1)*r2(1)
pv_bend(1,1) = pv_bend(1,1) + f3(1)*r3(1)
!p(1,2)
pvbend(1,2) = pvbend(1,2) + f1(1)*r1(2)
pvbend(1,2) = pvbend(1,2) + f2(1)*r2(2)
pvbend(1,2) = pvbend(1,2) + f3(1)*r3(2)
pv_bend(1,2) = pv_bend(1,2) + f1(1)*r1(2)
pv_bend(1,2) = pv_bend(1,2) + f2(1)*r2(2)
pv_bend(1,2) = pv_bend(1,2) + f3(1)*r3(2)
!p(1,3)
pvbend(1,3) = pvbend(1,3) + f1(1)*r1(3)
pvbend(1,3) = pvbend(1,3) + f2(1)*r2(3)
pvbend(1,3) = pvbend(1,3) + f3(1)*r3(3)
pv_bend(1,3) = pv_bend(1,3) + f1(1)*r1(3)
pv_bend(1,3) = pv_bend(1,3) + f2(1)*r2(3)
pv_bend(1,3) = pv_bend(1,3) + f3(1)*r3(3)
!p(2,1)
pvbend(2,1) = pvbend(2,1) + f1(2)*r1(1)
pvbend(2,1) = pvbend(2,1) + f2(2)*r2(1)
pvbend(2,1) = pvbend(2,1) + f3(2)*r3(1)
pv_bend(2,1) = pv_bend(2,1) + f1(2)*r1(1)
pv_bend(2,1) = pv_bend(2,1) + f2(2)*r2(1)
pv_bend(2,1) = pv_bend(2,1) + f3(2)*r3(1)
!p(2,2)
pvbend(2,2) = pvbend(2,2) + f1(2)*r1(2)
pvbend(2,2) = pvbend(2,2) + f2(2)*r2(2)
pvbend(2,2) = pvbend(2,2) + f3(2)*r3(2)
pv_bend(2,2) = pv_bend(2,2) + f1(2)*r1(2)
pv_bend(2,2) = pv_bend(2,2) + f2(2)*r2(2)
pv_bend(2,2) = pv_bend(2,2) + f3(2)*r3(2)
!p(2,3)
pvbend(2,3) = pvbend(2,3) + f1(2)*r1(3)
pvbend(2,3) = pvbend(2,3) + f2(2)*r2(3)
pvbend(2,3) = pvbend(2,3) + f3(2)*r3(3)
pv_bend(2,3) = pv_bend(2,3) + f1(2)*r1(3)
pv_bend(2,3) = pv_bend(2,3) + f2(2)*r2(3)
pv_bend(2,3) = pv_bend(2,3) + f3(2)*r3(3)
!p(3,1)
pvbend(3,1) = pvbend(3,1) + f1(3)*r1(1)
pvbend(3,1) = pvbend(3,1) + f2(3)*r2(1)
pvbend(3,1) = pvbend(3,1) + f3(3)*r3(1)
pv_bend(3,1) = pv_bend(3,1) + f1(3)*r1(1)
pv_bend(3,1) = pv_bend(3,1) + f2(3)*r2(1)
pv_bend(3,1) = pv_bend(3,1) + f3(3)*r3(1)
!p(3,2)
pvbend(3,2) = pvbend(3,2) + f1(3)*r1(2)
pvbend(3,2) = pvbend(3,2) + f2(3)*r2(2)
pvbend(3,2) = pvbend(3,2) + f3(3)*r3(2)
pv_bend(3,2) = pv_bend(3,2) + f1(3)*r1(2)
pv_bend(3,2) = pv_bend(3,2) + f2(3)*r2(2)
pv_bend(3,2) = pv_bend(3,2) + f3(3)*r3(2)
!p(3,3)
pvbend(3,3) = pvbend(3,3) + f1(3)*r1(3)
pvbend(3,3) = pvbend(3,3) + f2(3)*r2(3)
pvbend(3,3) = pvbend(3,3) + f3(3)*r3(3)
pv_bend(3,3) = pv_bend(3,3) + f1(3)*r1(3)
pv_bend(3,3) = pv_bend(3,3) + f2(3)*r2(3)
pv_bend(3,3) = pv_bend(3,3) + f3(3)*r3(3)
END SUBROUTINE get_pv_bend
!******************************************************************************
SUBROUTINE get_pv_torsion()
SUBROUTINE get_pv_torsion(f1,f2,f3,f4,r1,r2,r3,r4,pv_torsion)
IMPLICIT NONE
! Arguments
REAL ( dbl ), DIMENSION (3), INTENT ( IN ) :: f1, f2, f3, f4, r1, r2, r3, r4
REAL ( dbl ), DIMENSION ( :, : ), INTENT ( INOUT ) :: pv_torsion
!------------------------------------------------------------------------------
!p(1,1)
pv_torsion(1,1) = pv_torsion(1,1) + f1(1)*r1(1)
pv_torsion(1,1) = pv_torsion(1,1) + f2(1)*r2(1)
pv_torsion(1,1) = pv_torsion(1,1) + f3(1)*r3(1)
pv_torsion(1,1) = pv_torsion(1,1) + f4(1)*r4(1)
!p(1,2)
pv_torsion(1,2) = pv_torsion(1,2) + f1(1)*r1(2)
pv_torsion(1,2) = pv_torsion(1,2) + f2(1)*r2(2)
pv_torsion(1,2) = pv_torsion(1,2) + f3(1)*r3(2)
pv_torsion(1,2) = pv_torsion(1,2) + f4(1)*r4(2)
!p(1,3)
pv_torsion(1,3) = pv_torsion(1,3) + f1(1)*r1(3)
pv_torsion(1,3) = pv_torsion(1,3) + f2(1)*r2(3)
pv_torsion(1,3) = pv_torsion(1,3) + f3(1)*r3(3)
pv_torsion(1,3) = pv_torsion(1,3) + f4(1)*r4(3)
!p(2,1)
pv_torsion(2,1) = pv_torsion(2,1) + f1(2)*r1(1)
pv_torsion(2,1) = pv_torsion(2,1) + f2(2)*r2(1)
pv_torsion(2,1) = pv_torsion(2,1) + f3(2)*r3(1)
pv_torsion(2,1) = pv_torsion(2,1) + f4(2)*r4(1)
!p(2,2)
pv_torsion(2,2) = pv_torsion(2,2) + f1(2)*r1(2)
pv_torsion(2,2) = pv_torsion(2,2) + f2(2)*r2(2)
pv_torsion(2,2) = pv_torsion(2,2) + f3(2)*r3(2)
pv_torsion(2,2) = pv_torsion(2,2) + f4(2)*r4(2)
!p(2,3)
pv_torsion(2,3) = pv_torsion(2,3) + f1(2)*r1(3)
pv_torsion(2,3) = pv_torsion(2,3) + f2(2)*r2(3)
pv_torsion(2,3) = pv_torsion(2,3) + f3(2)*r3(3)
pv_torsion(2,3) = pv_torsion(2,3) + f4(2)*r4(3)
!p(3,1)
pv_torsion(3,1) = pv_torsion(3,1) + f1(3)*r1(1)
pv_torsion(3,1) = pv_torsion(3,1) + f2(3)*r2(1)
pv_torsion(3,1) = pv_torsion(3,1) + f3(3)*r3(1)
pv_torsion(3,1) = pv_torsion(3,1) + f4(3)*r4(1)
!p(3,2)
pv_torsion(3,2) = pv_torsion(3,2) + f1(3)*r1(2)
pv_torsion(3,2) = pv_torsion(3,2) + f2(3)*r2(2)
pv_torsion(3,2) = pv_torsion(3,2) + f3(3)*r3(2)
pv_torsion(3,2) = pv_torsion(3,2) + f4(3)*r4(2)
!p(3,3)
pv_torsion(3,3) = pv_torsion(3,3) + f1(3)*r1(3)
pv_torsion(3,3) = pv_torsion(3,3) + f2(3)*r2(3)
pv_torsion(3,3) = pv_torsion(3,3) + f3(3)*r3(3)
pv_torsion(3,3) = pv_torsion(3,3) + f4(3)*r4(3)
STOP 1990
END SUBROUTINE get_pv_torsion
!******************************************************************************

244
src/molecule_input.F
View file

@ -2,6 +2,20 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/molecule_input [1.0] *
!!
!! NAME
!! molecule_input
!!
!! FUNCTION
!!
!! AUTHOR
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE molecule_input
@ -21,29 +35,49 @@ MODULE molecule_input
PRIVATE
PUBLIC :: read_setup_section, read_molecule_section, charge, multiplicity
!!total charge of the system
! total charge of the system
REAL ( dbl ) :: charge
!!multiplicity of the system
! multiplicity of the system
INTEGER :: multiplicity
!!*****
!******************************************************************************
CONTAINS
!!>---------------------------------------------------------------------------!
!! SECTION: &setup ... &end !
!! !
!! charge charge !
!! multiplicity multiplicity !
!! moltype nmol_types !
!! molecule !
!! name molname !
!! nummol nmol !
!! nhcopt [molecule,global,massive,...] !
!! end molecule !
!! !
!!<---------------------------------------------------------------------------!
!******************************************************************************
!!****** molecule_input/read_setup_section [1.0] *
!!
!! NAME
!! read_setup_section
!!
!! FUNCTION
!!
!! AUTHOR
!!
!! MODIFICATION HISTORY
!! none
!!
!! NOTES
!!I---------------------------------------------------------------------------I
!!I SECTION: &setup ... &end I
!!I I
!!I charge charge I
!!I multiplicity multiplicity I
!!I moltype nmol_types I
!!I molecule I
!!I name molname I
!!I nummol nmol I
!!I nhcopt [molecule,global,massive,...] I
!!I end molecule I
!!I I
!!I---------------------------------------------------------------------------I
!!
!! SOURCE
!******************************************************************************
SUBROUTINE read_setup_section ( mol_setup, set_fn, globenv )
@ -154,24 +188,40 @@ SUBROUTINE read_setup_section ( mol_setup, set_fn, globenv )
END SUBROUTINE read_setup_section
!!*****
!******************************************************************************
!!****** cp2k/molecule_input [1.0] *
!!
!! NAME
!! molecule_input
!!
!! FUNCTION
!!
!! AUTHOR
!!
!! MODIFICATION HISTORY
!! none
!!
!! NOTES
!!I---------------------------------------------------------------------------I
!!I SECTION: &molecule ... &end I
!!I I
!!I moleculename I
!!I numat n I
!!I atom name label mass charge I
!!I bonds nbonds I
!!I bends nbends I
!!I torsion ntorsion I
!!I distcon ndcon I
!!I 3x3con n3x3con I
!!I basis name library nmax epsilon I
!!I drho name library nmax epsilon I
!!I end moleculename I
!!I I
!!I---------------------------------------------------------------------------I
!!
!! SOURCE
!******************************************************************************
!!>---------------------------------------------------------------------------!
!! SECTION: &molecule ... &end !
!! !
!! moleculename !
!! numat n !
!! atom name label mass charge !
!! bonds nbonds !
!! bends nbends !
!! torsion ntorsion !
!! distcon ndcon !
!! 3x3con n3x3con !
!! basis name library nmax epsilon !
!! drho name library nmax epsilon !
!! end moleculename !
!! !
!!<---------------------------------------------------------------------------!
SUBROUTINE read_molecule_section ( mol_setup, set_fn, globenv )
@ -275,23 +325,8 @@ SUBROUTINE read_molecule_section ( mol_setup, set_fn, globenv )
IF (ios /= 0 ) CALL stop_memory &
( 'read_molecule_section', 'bonds', 2 * na )
DO i = 1, na
ilen = 0
CALL cfield(at1,ilen)
i1 = str_search ( &
mol_setup ( imol ) % molpar % alabel, nat, at1 )
IF ( i1 == 0 ) &
CALL stop_parser ( 'read_molecule_section', &
'atom label in bonds not found' )
ilen = 0
CALL cfield(at2,ilen)
i2 = str_search ( &
mol_setup ( imol ) % molpar % alabel, nat, at2 )
IF ( i2 == 0 ) &
CALL stop_parser ( 'read_molecule_section', &
'atom label in bonds not found' )
mol_setup(imol) %molpar%bonds(1,i) = i1
mol_setup(imol) %molpar%bonds(2,i) = i2
mol_setup(imol) %molpar%bonds(1,i) = get_int()
mol_setup(imol) %molpar%bonds(2,i) = get_int()
END DO
END IF
@ -303,30 +338,9 @@ SUBROUTINE read_molecule_section ( mol_setup, set_fn, globenv )
IF (ios /= 0 ) CALL stop_memory &
( 'read_molecule_section', 'bends', 3 * na )
DO i = 1, na
ilen = 0
CALL cfield(at1,ilen)
i1 = str_search ( &
mol_setup(imol)%molpar%alabel,nat,at1)
IF ( i1 == 0 ) CALL stop_parser &
( 'read_molecule_section', &
'atom label in bends not found' )
ilen = 0
CALL cfield(at2,ilen)
i2 = str_search ( &
mol_setup(imol)%molpar%alabel,nat,at2)
IF ( i2 == 0 ) &
CALL stop_parser ( 'read_molecule_section', &
'atom label in bends not found' )
ilen = 0
CALL cfield(at3,ilen)
i3 = str_search ( &
mol_setup(imol)%molpar%alabel,nat,at3)
IF ( i3 == 0 ) &
CALL stop_parser ( 'read_molecule_section', &
'atom label in bends not found' )
mol_setup(imol) %molpar%bends(1,i) = i1
mol_setup(imol) %molpar%bends(2,i) = i2
mol_setup(imol) %molpar%bends(3,i) = i3
mol_setup(imol) %molpar%bends(1,i) = get_int()
mol_setup(imol) %molpar%bends(2,i) = get_int()
mol_setup(imol) %molpar%bends(3,i) = get_int()
END DO
END IF
@ -340,38 +354,10 @@ SUBROUTINE read_molecule_section ( mol_setup, set_fn, globenv )
IF ( ios /= 0 ) CALL stop_memory &
( 'read_molecule_section', 'torsions', 4 * na )
DO i = 1, na
ilen = 0
CALL cfield(at1,ilen)
i1 = str_search ( &
mol_setup(imol)%molpar%alabel,nat,at1)
IF ( i1 == 0 ) CALL stop_parser &
( 'read_molecule_section', &
'atom label in torsions not found' )
ilen = 0
CALL cfield(at2,ilen)
i2 = str_search ( &
mol_setup(imol)%molpar%alabel,nat,at2)
IF ( i2 == 0 ) CALL stop_parser &
( 'read_molecule_section', &
'atom label in torsions not found' )
ilen = 0
CALL cfield(at3,ilen)
i3 = str_search ( &
mol_setup(imol)%molpar%alabel,nat,at3)
IF ( i3 == 0 ) CALL stop_parser &
( 'read_molecule_section', &
'atom label in torsions not found' )
ilen = 0
CALL cfield(at4,ilen)
i4 = str_search ( &
mol_setup(imol)%molpar%alabel,nat,at4)
IF ( i4 == 0 ) CALL stop_parser &
( 'read_molecule_section', &
'atom label in torsions not found' )
mol_setup(imol) %molpar%torsions(1,i) = i1
mol_setup(imol) %molpar%torsions(2,i) = i2
mol_setup(imol) %molpar%torsions(3,i) = i3
mol_setup(imol) %molpar%torsions(4,i) = i4
mol_setup(imol) %molpar%torsions(1,i) = get_int()
mol_setup(imol) %molpar%torsions(2,i) = get_int()
mol_setup(imol) %molpar%torsions(3,i) = get_int()
mol_setup(imol) %molpar%torsions(4,i) = get_int()
END DO
END IF
@ -383,22 +369,8 @@ SUBROUTINE read_molecule_section ( mol_setup, set_fn, globenv )
IF ( ios /= 0 ) CALL stop_memory &
( 'read_molecule_section', 'ndc', 2 * na )
DO i = 1, na
ilen = 0
CALL cfield(at1,ilen)
i1 = str_search ( &
mol_setup(imol)%molpar%alabel,nat,at1)
IF ( i1 == 0 ) CALL stop_parser &
( 'read_molecule_section', &
'atom label in distcons not found' )
ilen = 0
CALL cfield(at2,ilen)
i2 = str_search ( &
mol_setup(imol)%molpar%alabel,nat,at2)
IF ( i2 == 0 ) CALL stop_parser &
( 'read_molecule_section', &
'atom label in distcons not found' )
mol_setup(imol) %molpar%ndc(1,i) = i1
mol_setup(imol) %molpar%ndc(2,i) = i2
mol_setup(imol) %molpar%ndc(1,i) = get_int()
mol_setup(imol) %molpar%ndc(2,i) = get_int()
END DO
END IF
@ -410,30 +382,9 @@ SUBROUTINE read_molecule_section ( mol_setup, set_fn, globenv )
IF (ios /= 0 ) CALL stop_memory &
( 'read_molecule_section', 'n3x3c', 3 * na )
DO i = 1, na
ilen = 0
CALL cfield(at1,ilen)
i1 = str_search ( &
mol_setup(imol)%molpar%alabel,nat,at1)
IF ( i1 == 0 ) CALL stop_parser &
( 'read_molecule_section', &
'atom label in g3x3con not found' )
ilen = 0
CALL cfield(at2,ilen)
i2 = str_search ( &
mol_setup(imol)%molpar%alabel,nat,at2)
IF ( i2 == 0 ) CALL stop_parser &
( 'read_molecule_section', &
'atom label in g3x3con not found' )
ilen = 0
CALL cfield(at3,ilen)
i3 = str_search ( &
mol_setup(imol)%molpar%alabel,nat,at3)
IF ( i3 == 0 ) CALL stop_parser &
( 'read_molecule_section', &
'atom label in g3x3con not found' )
mol_setup(imol) %molpar%n3x3c(1,i) = i1
mol_setup(imol) %molpar%n3x3c(2,i) = i2
mol_setup(imol) %molpar%n3x3c(3,i) = i3
mol_setup(imol) %molpar%n3x3c(1,i) = get_int()
mol_setup(imol) %molpar%n3x3c(2,i) = get_int()
mol_setup(imol) %molpar%n3x3c(3,i) = get_int()
END DO
END IF
@ -610,6 +561,9 @@ SUBROUTINE read_molecule_section ( mol_setup, set_fn, globenv )
END SUBROUTINE read_molecule_section
!!*****
!******************************************************************************
END MODULE molecule_input
!******************************************************************************

20
src/molecule_types.F
View file

@ -2,6 +2,21 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****s* cp2k/molecule_types [1.0] *
!!
!! NAME
!! molecule_types
!!
!! FUNCTION
!!
!! AUTHOR
!! CJM & JGH
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE molecule_types
@ -84,7 +99,7 @@ MODULE molecule_types
TYPE torsion_parameters_type
CHARACTER ( LEN = 20 ) :: type
REAL ( dbl ), DIMENSION (0:6) :: parm
REAL ( dbl ), DIMENSION (0:3) :: a
END TYPE torsion_parameters_type
TYPE linklist_torsions
@ -171,3 +186,6 @@ MODULE molecule_types
END TYPE molecule_structure_type
END MODULE molecule_types
!******************************************************************************
!!*****

44
src/nhc.F
View file

@ -2,6 +2,21 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/nhc [1.0] *
!!
!! NAME
!! nhc
!!
!! FUNCTION
!!
!! AUTHOR
!! CJM
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE nhc
@ -21,11 +36,28 @@ MODULE nhc
REAL ( dbl ), TARGET :: scale
!!*****
!******************************************************************************
CONTAINS
!******************************************************************************
! creates the thermostatting maps use global_mod
!!****** nhc/nhc_control [1.0] *
!!
!! NAME
!! nhc_control
!!
!! FUNCTION
!! Creates the thermostatting maps
!!
!! AUTHOR
!! CJM
!!
!! MODIFICATION HISTORY
!! 29-Nov-00 (JGH) correct counting of DOF if constraints are off
!!
!! SOURCE
!******************************************************************************
SUBROUTINE nhc_control ( simpar, molecule, mol_setup, nhcp, globenv )
@ -201,6 +233,11 @@ SUBROUTINE nhc_control ( simpar, molecule, mol_setup, nhcp, globenv )
globenv%mepos)
END DO
END IF
IF ( .NOT. simpar % constraint ) THEN
const_mol = 0._dbl
tot_const = 0._dbl
END IF
! count the number of atoms on this processor
natoms = sum(molecule ( : )%natoms_mol)
@ -620,6 +657,9 @@ SUBROUTINE nhc_control ( simpar, molecule, mol_setup, nhcp, globenv )
END SUBROUTINE nhc_control
!!*****
!******************************************************************************
END MODULE nhc
!******************************************************************************

156
src/nose.F
View file

@ -2,6 +2,21 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/nose [1.0] *
!!
!! NAME
!! nose
!!
!! FUNCTION
!!
!! AUTHOR
!! CJM
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE nose
@ -59,8 +74,28 @@ MODULE nose
REAL ( dbl ), POINTER :: s_kin ( : )
END TYPE extended_parameters_type
!!*****
!******************************************************************************
CONTAINS
!******************************************************************************
!!****** nose/lnhc [1.0] *
!!
!! NAME
!! lnhc
!!
!! FUNCTION
!!
!! AUTHOR
!!
!! CREATION DATE
!! 14-NOV-2000
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
SUBROUTINE lnhc ( nhcp, pnode, group )
@ -194,6 +229,24 @@ SUBROUTINE lnhc ( nhcp, pnode, group )
END SUBROUTINE lnhc
!!*****
!******************************************************************************
!!****** nose/lnhcp [1.0] *
!!
!! NAME
!! lnhcp
!!
!! FUNCTION
!!
!! AUTHOR
!!
!! CREATION DATE
!! 14-NOV-2000
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
SUBROUTINE lnhcp ( nhcp, molecule, pnode, pv_virial, pv_kin, pv_const, pext, &
@ -392,6 +445,24 @@ SUBROUTINE lnhcp ( nhcp, molecule, pnode, pv_virial, pv_kin, pv_const, pext, &
END SUBROUTINE lnhcp
!!*****
!******************************************************************************
!!****** nose/lnhcpf [1.0] *
!!
!! NAME
!! lnhcpf
!!
!! FUNCTION
!!
!! AUTHOR
!!
!! CREATION DATE
!! 14-NOV-2000
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
SUBROUTINE lnhcpf ( nhcp, molecule, pnode, pv_virial, pv_kin, pv_const, pext, &
@ -643,10 +714,29 @@ SUBROUTINE lnhcpf ( nhcp, molecule, pnode, pv_virial, pv_kin, pv_const, pext, &
flops = flops*1.E-6_dbl
CALL timestop(flops,handle)
END SUBROUTINE lnhcpf
!!*****
!******************************************************************************
!!****** nose/initial_nhc [1.0] *
!!
!! NAME
!! initial_nhc
!!
!! FUNCTION
!! Initializes the NHC velocities to the Maxwellian distribution
!!
!! AUTHOR
!!
!! CREATION DATE
!! 14-NOV-2000
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
SUBROUTINE initial_nhc(nhcp,tau_nhc,temp_ext, globenv )
! Initializes the NHC velocities to the Maxwellian distribution
IMPLICIT NONE
TYPE (extended_parameters_type ), INTENT ( INOUT ) :: nhcp
TYPE ( global_environment_type ), INTENT ( INOUT ) :: globenv
@ -680,15 +770,12 @@ SUBROUTINE initial_nhc(nhcp,tau_nhc,temp_ext, globenv )
IF (nhcp % dis_type=='DIS_REP') CALL mp_sum(number,globenv % group)
#endif
temp = 2.0_dbl*akin/REAL(number)
! scale velocities to get the correct initial temperature
temp = 2.0_dbl*akin/REAL(number)
temp = sqrt(temp_ext/temp)
DO i = 1, nhcp % num_nhc
DO j = 1, nhcp % nhc_len
nhcp % nvt(j,i) % v = sqrt(temp_ext/temp)*nhcp % nvt(j,i) % v
! dbg
! nhcp % nvt(j,i) % v = 0._dbl
! dbg
nhcp % nvt(j,i) % v = temp * nhcp % nvt(j,i) % v
nhcp % nvt(j,i) % eta = 0.0_dbl
END DO
END DO
@ -704,9 +791,26 @@ SUBROUTINE initial_nhc(nhcp,tau_nhc,temp_ext, globenv )
END SUBROUTINE initial_nhc
!!*****
!******************************************************************************
!!****** nose/initial_npt [1.0] *
!!
!! NAME
!! initial_npt
!!
!! FUNCTION
!! Initializes the NHC velocities to the Maxwellian distribution
!!
!! AUTHOR
!!
!! CREATION DATE
!! 14-NOV-2000
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
! Initializes the NHC velocities to the Maxwellian distribution
SUBROUTINE initial_npt ( nhcp, tau_cell, temp_ext, deth, nfree, ensemble, &
globenv )
@ -763,15 +867,12 @@ SUBROUTINE initial_npt ( nhcp, tau_cell, temp_ext, deth, nfree, ensemble, &
END DO
number = SIZE ( nhcp % npt, 1 ) * SIZE ( nhcp % npt, 2 )
temp = 2.0_dbl * akin / REAL ( number, dbl )
! scale velocities to get the correct initial temperature
temp = 2.0_dbl * akin / REAL ( number, dbl )
temp = sqrt(temp_ext/temp)
DO i = 1, size(nhcp % npt,1)
DO j = 1, size(nhcp % npt,2)
nhcp % npt(j,i) % v = sqrt(temp_ext/temp)*nhcp % npt(j,i) % v
! dbg
! nhcp % npt(j,i) % v = 0._dbl
! dbg
nhcp % npt(j,i) % v = temp * nhcp % npt(j,i) % v
END DO
END DO
@ -780,10 +881,27 @@ SUBROUTINE initial_npt ( nhcp, tau_cell, temp_ext, deth, nfree, ensemble, &
END SUBROUTINE initial_npt
!!*****
!******************************************************************************
!!****** nose/yoshida_coef [1.0] *
!!
!! NAME
!! yoshida_coef
!!
!! FUNCTION
!! This lists the coefficients for the Yoshida method (higher
!! order integrator used in NVT)
!!
!! AUTHOR
!!
!! CREATION DATE
!! 14-NOV-2000
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
! this lists the coefficients for the Yoshida method (higher
! order integrator used in NVT)
SUBROUTINE yoshida_coef ( nhcp, dt )
@ -854,7 +972,9 @@ SUBROUTINE yoshida_coef ( nhcp, dt )
END SUBROUTINE yoshida_coef
!!*****
!******************************************************************************
END MODULE nose
!******************************************************************************

15
src/nrutil.F
View file

@ -1,3 +1,17 @@
!!****** cp2k/nrutil [1.0] *
!!
!! NAME
!! nrutil
!!
!! FUNCTION
!!
!! AUTHOR
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE nrutil
USE kinds, ONLY : dbl
IMPLICIT NONE
@ -73,6 +87,7 @@
MODULE PROCEDURE put_diag_rv, put_diag_r
END INTERFACE
!!*****
CONTAINS
!BL
SUBROUTINE array_copy_r(src,dest,n_copied,n_not_copied)

105
src/output_utilities.F
View file

@ -2,10 +2,22 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
! Purpose: Output utilities.
! History: - Creation (06.10.1999, Matthias Krack)
! - CP2K by JH 21.08.2000
!!****** cp2k/output_utilities [1.0] *
!!
!! NAME
!! output_utilities
!!
!! AUTHOR
!! Matthias Krack
!!
!! CREATION DATE
!! 06.10.1999
!!
!! MODIFICATION HISTORY
!! CP2K Version by JGH 21.08.2000
!!
!! SOURCE
!******************************************************************************
MODULE output_utilities
@ -16,27 +28,40 @@ MODULE output_utilities
INTEGER, PARAMETER :: out = 6
!!*****
!******************************************************************************
CONTAINS
!******************************************************************************
! Purpose: Write a message string to the logical output unit lunit.
! History: - Creation (28.08.1996, Matthias Krack)
! - Translated to Fortran 90/95 (07.10.1999, Matthias Krack)
! - CP2K by JH 21.08.2000
!!****** output_utilities/print_message [1.0] *
!!
!! NAME
!! print_message
!!
!! AUTHOR
!! Matthias Krack
!!
!! CREATION DATE
!! 28.08.1996
!!
!! MODIFICATION HISTORY
!! - Creation (28.08.1996, Matthias Krack)
!! - Translated to Fortran 90/95 (07.10.1999, Matthias Krack)
!! - CP2K by JH 21.08.2000
!!
!! NOTES
!! after : Number of empty lines after the message.
!! before : Number of empty lines before the message.
!! declev : Decoration level (0,1,2, ... star lines).
!! lunit : Number of the logical output unit.
!! message: String with the message text.
!!
!! SOURCE
!******************************************************************************
SUBROUTINE print_message ( message, lunit, declev, before, after )
! ***************************************************************************
! after : Number of empty lines after the message.
! before : Number of empty lines before the message.
! declev : Decoration level (0,1,2, ... star lines).
! lunit : Number of the logical output unit.
! message: String with the message text.
! ***************************************************************************
! Arguments
CHARACTER ( LEN = * ), INTENT ( IN ) :: message
INTEGER, INTENT ( IN ), OPTIONAL :: after, before, declev, lunit
@ -50,8 +75,6 @@ SUBROUTINE print_message ( message, lunit, declev, before, after )
!------------------------------------------------------------------------------
! IF (my_pe /= io_pe) RETURN
IF ( PRESENT ( after ) ) THEN
blank_lines_after = MAX ( after, 0 )
ELSE
@ -180,23 +203,34 @@ SUBROUTINE print_message ( message, lunit, declev, before, after )
END SUBROUTINE print_message
!!*****
!******************************************************************************
!!****** output_utilities/print_warning [1.0] *
!!
!! NAME
!! print_warning
!!
!! AUTHOR
!! Matthias Krack
!!
!! CREATION DATE
!! 28.08.1996
!!
!! MODIFICATION HISTORY
!! - Creation (28.08.1996, Matthias Krack)
!! - Translated to Fortran 90/95 (07.10.1999, Matthias Krack)
!! - CP2K by JH 21.08.2000
!!
!! NOTES
!! message: Message text.
!! routine: Name of the calling routine.
!! lunit : Output unit
!!
!! SOURCE
!******************************************************************************
! Purpose: Print a warning message.
! History: - Creation (28.08.1996, Matthias Krack)
! - Translated to Fortran 90/95 (07.10.1999, Matthias Krack)
! - CP2K by JH 21.08.2000
SUBROUTINE print_warning ( routine, message, lunit )
! ***************************************************************************
! message: Message text.
! routine: Name of the calling routine.
! lunit : Output unit
! ***************************************************************************
! Arguments
CHARACTER ( LEN = * ), INTENT ( IN ) :: message, routine
INTEGER, INTENT ( IN ), OPTIONAL :: lunit
@ -215,6 +249,9 @@ SUBROUTINE print_warning ( routine, message, lunit )
END SUBROUTINE print_warning
!!*****
!******************************************************************************
END MODULE output_utilities
!******************************************************************************

18
src/pair_potential.F
View file

@ -2,6 +2,21 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/pair_potential [1.0] *
!!
!! NAME
!! pair_potential
!!
!! FUNCTION
!!
!! AUTHOR
!! CJM
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE pair_potential
@ -55,6 +70,9 @@ MODULE pair_potential
! splines for interaction n
TYPE ( spline_data ), DIMENSION ( : ), POINTER, SAVE :: spline_table
!!*****
!******************************************************************************
CONTAINS
!******************************************************************************

20
src/parallel.F
View file

@ -2,6 +2,21 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/parallel [1.0] *
!!
!! NAME
!! parallel
!!
!! FUNCTION
!!
!! AUTHOR
!! JGH
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE parallel
@ -12,6 +27,8 @@ MODULE parallel
PRIVATE
PUBLIC :: start_parallel, end_parallel
!******************************************************************************
CONTAINS
!******************************************************************************
@ -47,3 +64,6 @@ END SUBROUTINE end_parallel
!******************************************************************************
END MODULE parallel
!******************************************************************************
!!*****

71
src/parser.F
View file

@ -2,30 +2,48 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!-----------------------------------------------------------------------------!
!!> F R E E - F O R M A T I N P U T R O U T I N E S !
!! !
!! Function INSCAN(iunit,label) : looks for an input dection started by !
!! label un file associated with unit iunit !
!! Subroutine READ_LINE(iunit) : reads an input line, physical end of lines !
!! can be masked by a & character, no further !
!! characters are alowed to the right !
!! Function TEST_NEXT() : tests contains of next input field, return values!
!! are : E (blank); C (text); N (number); !
!! K (bracket); S (! comment); X (&END); !
!! O all other !
!! Subroutine CFIELD(string,length) : extracts the next field from the !
!! input line, at most length charcters !
!! are returned in string, if length is 0 !
!! the full field is returned, possible !
!! field seperators are: blank, comma, !
!! colon, semicolon, bracket, equal sign !
!! Subroutine P_ERROR : writes out current input line and current position !
!! within the line, does not stop code !
!! Function GET_INT() : read an integer from the input line !
!! Function GET_REAL() : read an real value from the input line !
!!<Function GET_ARRAY(I,J) : read an array field [I:J] or [I] !
!-----------------------------------------------------------------------------!
!!****** cp2k/parser [1.0] *
!!
!! NAME
!! parser
!!
!! FUNCTION
!! Provides all the functions needed to parse an input section
!!
!! AUTHOR
!! JGH
!!
!! MODIFICATION HISTORY
!! none
!!
!! NOTES
!!----------------------------------------------------------------------------
!! F R E E - F O R M A T I N P U T R O U T I N E S
!!
!! Function INSCAN(iunit,label) : looks for an input dection started by
!! label un file associated with unit iunit
!! Subroutine READ_LINE(iunit) : reads an input line, physical end of lines
!! can be masked by a & character, no further
!! characters are alowed to the right
!! Function TEST_NEXT() : tests contains of next input field, return values
!! are : E (blank); C (text); N (number);
!! K (bracket); S (! comment); X (&END);
!! O all other
!! Subroutine CFIELD(string,length) : extracts the next field from the
!! input line, at most length charcters
!! are returned in string, if length is 0
!! the full field is returned, possible
!! field seperators are: blank, comma,
!! colon, semicolon, bracket, equal sign
!! Subroutine P_ERROR : writes out current input line and current position
!! within the line, does not stop code
!! Function GET_INT() : read an integer from the input line
!! Function GET_REAL() : read an real value from the input line
!! Function GET_ARRAY(I,J) : read an array field [I:J] or [I]
!!----------------------------------------------------------------------------
!!
!! SOURCE
!-----------------------------------------------------------------------------
MODULE parser
@ -52,6 +70,8 @@ MODULE parser
TYPE(global_environment_type) :: parenv
LOGICAL :: parallel_io, ionode
!******************************************************************************
CONTAINS
!******************************************************************************
@ -785,3 +805,6 @@ END SUBROUTINE stop_parser
!******************************************************************************
END MODULE parser
!******************************************************************************
!!*****

18
src/particle_lists.F
View file

@ -2,6 +2,21 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/particle_lists [1.0] *
!!
!! NAME
!! particle_lists
!!
!! FUNCTION
!!
!! AUTHOR
!! CJM
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE particle_lists
@ -21,6 +36,9 @@ MODULE particle_lists
TYPE ( particle_list_type ), POINTER :: next
END TYPE particle_list_type
!!*****
!******************************************************************************
CONTAINS
!******************************************************************************

18
src/particle_types.F
View file

@ -2,6 +2,21 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****s* cp2k/particle_types [1.0] *
!!
!! NAME
!! particle_types
!!
!! FUNCTION
!!
!! AUTHOR
!! CJM
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE particle_types
@ -29,3 +44,6 @@ MODULE particle_types
END TYPE particle_type
END MODULE particle_types
!******************************************************************************
!!*****

40
src/periodic_table.F
View file

@ -2,8 +2,22 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!! Periodic Table related data definitions
!!****** cp2k/periodic_table [1.0] *
!!
!! NAME
!! periodic_table
!!
!! FUNCTION
!! Periodic Table related data definitions
!!
!! AUTHOR
!! JGH
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE periodic_table
@ -30,8 +44,27 @@ MODULE periodic_table
INTEGER, PARAMETER :: nelem = 106
TYPE ( atom ) :: ptable ( 0:nelem )
!!*****
!******************************************************************************
CONTAINS
!******************************************************************************
!!****** periodic_table/init_periodic_table [1.0] *
!!
!! NAME
!! init_periodic_table
!!
!! FUNCTION
!! Initialization of Periodic Table related data
!!
!! AUTHOR
!! JGH
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
SUBROUTINE init_periodic_table()
@ -1227,6 +1260,9 @@ SUBROUTINE init_periodic_table()
END SUBROUTINE init_periodic_table
!!*****
!******************************************************************************
END MODULE periodic_table
!******************************************************************************

105
src/physcon.F
View file

@ -2,33 +2,48 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/physcon [1.0] *
!!
!! NAME
!! physcon
!!
!! FUNCTION
!! Definition of physical constants.
!!
!! List of variables:
!!
!! a_bohr : Bohr radius [m].
!! a_fine : Fine-structure constant.
!! angstrom : Conversion factor [Bohr] -> [Angstrom]
!! bohr : Conversion factor [Angstrom] -> [Bohr].
!! c_light : Speed of light in vacuum [m/s].
!! e_charge : Elementary charge [C].
!! e_mass : Electron mass [kg].
!! a_mass : Atomic mass unit[kg].
!! esu : Conversion factors [Hartree] -> [esu].
!! evolt : Conversion factor [Hartree] -> [eV].
!! h_planck : Planck constant [J*s].
!! hz : Conversion factor [Hartree] -> [Hz].
!! kcalmol : Conversion factor [Hartree] -> [kcal/mol].
!! kjmol : Conversion factor [Hartree] -> [kJ/mol].
!! joule : Conversion factor [Hartree] -> [J].
!! mhz : Conversion factor [Hartree] -> [MHz].
!! mu_perm : Permeability of vacuum [N/A**2].
!! boltzmann : Boltzman constant [J/Kelvin].
!! n_avogadro: Avogrado constant [1/mol].
!! rydberg : Rydberg constant [1/m].
!! vibfac : Conversion factor [Hartree/Bohr**2] -> [1/cm].
!! wavenum : Conversion factor [Hartree] -> [1/cm] (wave numbers).
!!
!! AUTHOR
!! Matthias Krack
!!
!! MODIFICATION HISTORY
!! Adapted for use in CP2K (JGH)
!!
!! SOURCE
!******************************************************************************
!
!!>Definition of physical constants.
!!
!! List of variables:
!!
!! a_bohr : Bohr radius [m].
!! a_fine : Fine-structure constant.
!! angstrom : Conversion factor [Bohr] -> [Angstrom]
!! bohr : Conversion factor [Angstrom] -> [Bohr].
!! c_light : Speed of light in vacuum [m/s].
!! e_charge : Elementary charge [C].
!! e_mass : Electron mass [kg].
!! a_mass : Atomic mass unit[kg].
!! esu : Conversion factors [Hartree] -> [esu].
!! evolt : Conversion factor [Hartree] -> [eV].
!! h_planck : Planck constant [J*s].
!! hz : Conversion factor [Hartree] -> [Hz].
!! kcalmol : Conversion factor [Hartree] -> [kcal/mol].
!! kjmol : Conversion factor [Hartree] -> [kJ/mol].
!! joule : Conversion factor [Hartree] -> [J].
!! mhz : Conversion factor [Hartree] -> [MHz].
!! mu_perm : Permeability of vacuum [N/A**2].
!! boltzmann : Boltzman constant [J/Kelvin].
!! n_avogadro: Avogrado constant [1/mol].
!! rydberg : Rydberg constant [1/m].
!! vibfac : Conversion factor [Hartree/Bohr**2] -> [1/cm].
!!<wavenum : Conversion factor [Hartree] -> [1/cm] (wave numbers).
MODULE physcon
@ -50,12 +65,36 @@ MODULE physcon
REAL ( dbl ), DIMENSION (maxmom) :: esu
!!*****
!******************************************************************************
CONTAINS
!******************************************************************************
!!****** physcon/init_physcon [1.0] *
!!
!! NAME
!! init_physcon
!!
!! FUNCTION
!!
!! AUTHOR
!!
!! CREATION DATE
!! 14-NOV-2000
!!
!! MODIFICATION HISTORY
!! none
!!
!! NOTES
!! Fundamental physical constants (SI units)
!! Literature: The 1986 CODATA Recommended Values of the
!! Fundamental Physical Constants,
!! E. R. Cohen, B. N. Taylor,
!! J. Res. Natl. Bur. Stand. 92(2), 1 (1987)
!! SOURCE
!******************************************************************************
!! Initialization of physical constants.
SUBROUTINE init_physcon()
USE mathconstants, ONLY : pi
@ -67,13 +106,6 @@ SUBROUTINE init_physcon()
INTEGER :: i
!-----------------------------------------------------------------------------!
!
! *** Fundamental physical constants (SI units) ***
!
! *** Literature: The 1986 CODATA Recommended Values of the ***
! *** Fundamental Physical Constants, ***
! *** E. R. Cohen, B. N. Taylor, ***
! *** J. Res. Natl. Bur. Stand. 92(2), 1 (1987) ***
! *** Speed of light in vacuum [m/s] ***
c_light = 2.99792458E+8_dbl
@ -216,6 +248,9 @@ SUBROUTINE print_physcon ( iunit )
END SUBROUTINE print_physcon
!!*****
!******************************************************************************
END MODULE physcon
!******************************************************************************

41
src/pme.F
View file

@ -2,6 +2,20 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/pme [1.0] *
!!
!! NAME
!! pme
!!
!! FUNCTION
!!
!! AUTHOR
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE pme
@ -31,6 +45,9 @@ MODULE pme
PRIVATE
PUBLIC :: pme_evaluate, pme_setup
!!*****
!******************************************************************************
CONTAINS
!******************************************************************************
@ -182,17 +199,19 @@ SUBROUTINE pme_evaluate ( dg, fg_coulomb, vg_coulomb, pv_g, box, &
ELSE
size_of_rhos_r_save = SIZE ( rhos_r_save_dbl ) * bytes_per_real
END IF
IF ( size_of_rhos_r_save /= previous_size_of_rhos_r_save ) THEN
IF ( size_of_rhos_r_save / 1024**2 > 4 ) THEN
WRITE ( 6, '( /, A, I4, A, / )' ) " DG| SIZE(rho) = ", &
size_of_rhos_r_save / 1024**2, " MB"
ELSE
WRITE ( 6, '( /, A, I4, A, / )' ) " DG| SIZE(rho) = ", &
size_of_rhos_r_save / 1024, " kB"
END IF
previous_size_of_rhos_r_save = size_of_rhos_r_save
END IF
!
! this code has to be changed in order to be compliant with the output
! unit set and with parallel computing
! IF ( size_of_rhos_r_save /= previous_size_of_rhos_r_save ) THEN
! IF ( size_of_rhos_r_save / 1024**2 > 4 ) THEN
! WRITE ( 6, '( /, A, I4, A, / )' ) " DG| SIZE(rho) = ", &
! size_of_rhos_r_save / 1024**2, " MB"
! ELSE
! WRITE ( 6, '( /, A, I4, A, / )' ) " DG| SIZE(rho) = ", &
! size_of_rhos_r_save / 1024, " kB"
! END IF
! previous_size_of_rhos_r_save = size_of_rhos_r_save
! END IF
lb_pbc ( : ) = bounds_b ( 1, : ) + bounds_s ( 1, : )
ub_pbc ( : ) = bounds_b ( 2, : ) + bounds_s ( 2, : ) + 2

17
src/pw_grid_types.F
View file

@ -2,6 +2,20 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****s* cp2k/pw_grid_types [1.0] *
!!
!! NAME
!! pw_grid_types
!!
!! FUNCTION
!!
!! AUTHOR
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE pw_grid_types
@ -49,3 +63,6 @@ MODULE pw_grid_types
INTEGER, PARAMETER :: HALFSPACE = 211, FULLSPACE = 212
END MODULE pw_grid_types
!******************************************************************************
!!*****

24
src/pw_grids.F
View file

@ -2,10 +2,23 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!! This module defines the grid data type and some basic operations on it
!!
!!@Author apsi
!!@Author Christopher Mundy
!!****** cp2k/pw_grids [1.0] *
!!
!! NAME
!! pw_grids
!!
!! FUNCTION
!! This module defines the grid data type and some basic operations on it
!!
!! AUTHOR
!! apsi
!! Christopher Mundy
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE pw_grids
@ -21,6 +34,9 @@ MODULE pw_grids
PRIVATE
PUBLIC :: pw_grid_setup, pw_grid_change, pw_find_cutoff
!!*****
!******************************************************************************
CONTAINS
!******************************************************************************

17
src/pw_types.F
View file

@ -2,6 +2,20 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****s* cp2k/pw_types [1.0] *
!!
!! NAME
!! pw_types
!!
!! FUNCTION
!!
!! AUTHOR
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE pw_types
@ -39,6 +53,9 @@ MODULE pw_types
!! Flags for the structure member 'in_space'
INTEGER, PARAMETER :: NOSPACE = 371, REALSPACE = 372, RECIPROCALSPACE = 373
!!*****
!******************************************************************************
CONTAINS
!******************************************************************************

17
src/pws.F
View file

@ -2,6 +2,20 @@
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/pws [1.0] *
!!
!! NAME
!! pws
!!
!! FUNCTION
!!
!! AUTHOR
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE pws
@ -29,6 +43,9 @@ MODULE pws
fft_wrap_pw1c2, fft_wrap_pw1pw2
END INTERFACE
!!*****
!******************************************************************************
CONTAINS
!******************************************************************************

204
src/radial_grids.F Normal file
View file

@ -0,0 +1,204 @@
!-----------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!****** cp2k/radial_grids [1.0] *
!!
!! NAME
!! radial_grids
!!
!! FUNCTION
!! Radial grid transformation from
!! D. Andrae and J. Hinze, IJQC Vol. 63 pp 65 (1997)
!!
!! AUTHOR
!! JGH 17-NOV-2000
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE radial_grids
USE kinds, ONLY : dbl
USE string_utilities, ONLY : uppercase
USE stop_program, ONLY : stop_prg, stop_memory
IMPLICIT NONE
PRIVATE
PUBLIC :: radial_grid_type, init_radial_grid
TYPE radial_grid_type
CHARACTER (LEN=30) :: type
CHARACTER (LEN=60) :: form
INTEGER :: n
REAL(dbl) :: h
LOGICAL :: closed
REAL(dbl) :: z,b(5)
REAL(dbl) :: rmax
REAL(dbl) :: q
REAL(dbl), DIMENSION ( : ), POINTER :: r, w, wdd
END TYPE radial_grid_type
!!*****
!-----------------------------------------------------------------------------!
CONTAINS
!-----------------------------------------------------------------------------!
!!****** radial_grids/init_radial_grid [1.0] *
!!
!! NAME
!! init_radial_grid
!!
!! FUNCTION
!! Initialization routine for radial grids
!!
!! AUTHOR
!! JGH 17-NOV-2000
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!-----------------------------------------------------------------------------!
SUBROUTINE init_radial_grid ( rg, type, n, z, b, rmax )
!Arguments
TYPE ( radial_grid_type ), INTENT ( OUT ) :: rg
CHARACTER ( LEN = * ), INTENT ( IN ) :: type
INTEGER, INTENT ( IN ) :: n
REAL ( dbl ), INTENT ( IN ) :: z
REAL ( dbl ), DIMENSION ( : ), INTENT ( IN ), OPTIONAL :: b
REAL ( dbl ), INTENT ( IN ), OPTIONAL :: rmax
!Locals
REAL ( dbl ) :: ab, t, tm, t0, s, p, abt, ts2, tt
INTEGER :: ierr, m, i
!------------------------------------------------------------------------------
ALLOCATE ( rg % r ( 0:n+1 ), STAT = ierr )
IF ( ierr /= 0 ) call stop_memory ( "init_radial_grid", "rg % r", n+2 )
ALLOCATE ( rg % w ( 0:n+1 ), STAT = ierr )
IF ( ierr /= 0 ) call stop_memory ( "init_radial_grid", "rg % w", n+2 )
ALLOCATE ( rg % wdd ( 0:n+1 ), STAT = ierr )
IF ( ierr /= 0 ) call stop_memory ( "init_radial_grid", "rg % wdd", n+2 )
rg % type = type
CALL uppercase ( rg % type )
rg % n = n
rg % z = z
rg % h = 1._dbl / REAL ( n + 1, dbl )
IF ( PRESENT ( rmax ) ) THEN
rg % rmax = rmax
ELSE
rg % rmax = 100._dbl
END IF
IF ( PRESENT ( b ) ) THEN
rg % b = 0._dbl
m = SIZE ( b )
rg % b ( 1:m ) = b ( 1:m )
ELSE
rg % b = 1._dbl
END IF
SELECT CASE ( rg % type )
CASE DEFAULT
call stop_prg ( "init_radial_grid", "Variable transformation unknown" )
CASE ( "RATIONAL" )
rg % form = "(b/Z)[Ts/(1-Ts)]"
t0 = 0._dbl
tm = rg % z * rg % rmax / ( rg % z * rg % rmax + rg % b(1) )
t = tm - t0
ab = rg % z / rg % b(1)
abt = SQRT ( ab / t )
DO i = 0, n+1
s = REAL ( i, dbl ) / REAL ( n + 1, dbl )
rg % r ( i ) = ( t * s ) / ( 1._dbl - t * s ) / ab
rg % w ( i ) = abt * ( 1._dbl - t * s )
END DO
rg % wdd ( 0:n+1 ) = 0._dbl
rg % q = 1._dbl
rg % closed = .TRUE.
CASE ( "ALGEBRAIC SQUARE ROOT" )
rg % form = "(b/Z)[Ts/SQRT(1-(Ts)^2)]"
ab = rg % z / rg % b(1)
t0 = 0._dbl
tm = SIN ( ATAN ( ab * rg % rmax ) )
t = tm - t0
abt = SQRT ( ab / t )
DO i = 0, n+1
s = REAL ( i, dbl ) / REAL ( n + 1, dbl )
ts2 = t * t * s * s
rg % r ( i ) = ( t * s ) / SQRT ( 1._dbl - ts2 ) / ab
rg % w ( i ) = abt * ( 1._dbl - ts2 )**(0.75_dbl)
rg % wdd ( i ) = -0.75_dbl * t * t * abt * &
( 2._dbl - ts2 ) / ( 1._dbl - ts2 )**(1.25_dbl)
END DO
rg % q = 1._dbl
rg % closed = .TRUE.
CASE ( "LINEAR" )
rg % form = "Rmax*s"
DO i = 0, n+1
s = REAL ( i, dbl ) / REAL ( n + 1, dbl )
rg % r ( i ) = s * rmax
END DO
rg % w ( 0:n+1 ) = SQRT ( 1._dbl / rg % rmax )
rg % wdd ( 0:n+1 ) = 0._dbl
rg % q = 1._dbl
rg % closed = .FALSE.
CASE ( "POWER" )
rg % form = "(b/Z)[(Ts/b^p + 1)^(1/p) - 1]"
p = rg % b(2)
ab = rg % z / rg % b(1)
t0 = rg % b(1) ** p
tm = ( z * rg % rmax + rg % b(1) ) ** p
t = tm - t0
abt = SQRT ( p * rg % b(1)**p * ab / t )
tt = ( ( t * s ) / rg % b(1) ** p + 1._dbl )
DO i = 0, n+1
s = REAL ( i, dbl ) / REAL ( n + 1, dbl )
rg % r ( i ) = ( tt ** (1._dbl/p) - 1._dbl ) / ab
rg % w ( i ) = abt * tt ** ((p-1._dbl)/(2._dbl*p))
rg % wdd ( i ) = -0.25_dbl * abt * ( p - 1._dbl ) * ( p + 1._dbl ) * &
tt ** ((p-1._dbl)/(2._dbl*p)-2._dbl)
END DO
rg % q = 1._dbl
rg % closed = .FALSE.
CASE ( "LOG" )
rg % form = "(b/Z)[EXP(Ts) - 1]"
ab = rg % z / rg % b(1)
t0 = log ( rg % b(1) )
tm = log ( z * rg % rmax + rg % b(1) )
t = tm - t0
abt = SQRT ( ab / t )
DO i = 0, n+1
s = REAL ( i, dbl ) / REAL ( n + 1, dbl )
rg % r ( i ) = ( EXP ( t * s ) - 1._dbl ) / ab
rg % w ( i ) = abt * EXP ( -0.5_dbl * t * s )
rg % wdd ( i ) = 0.25_dbl * t * t * abt * EXP ( -0.5_dbl * t * s )
END DO
rg % q = 1._dbl
rg % closed = .FALSE.
END SELECT
END SUBROUTINE init_radial_grid
!!*****
!******************************************************************************
END MODULE radial_grids
!******************************************************************************

Some files were not shown because too many files have changed in this diff Show more