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

update to version 01/12/2000

Browse source 
svn-origin-rev: 9
This commit is contained in:
Jürg Hutter 2001-06-26 17:51:54 +00:00
parent 9a71617f7f
commit e221094355
46 changed files with 9341 additions and 1524 deletions
Download patch file Download diff file Expand all files Collapse all files

1
src/MACHINEDEFS.DEC
View file

@ -6,6 +6,7 @@ CPP = cpp
FC = f95 -free
FCfixed = f95 -fixed
LDR = f95
AR = ar -r
CDEFS = -D__DEC -D__FFTW

1
src/MACHINEDEFS.IBM
View file

@ -6,6 +6,7 @@ CPP = cpp
FC = xlf90
FCfixed = xlf90 -qfixed
LDR = xlf90
AR = ar -r
CDEFS = -D__AIX -D__FFTW

1
src/MACHINEDEFS.PGI
View file

@ -6,6 +6,7 @@ CPP = cpp
FC = pgf90 -Mfree
FCfixed = pgf90 -Mfixed
LDR = pgf90
AR = ar -r
CDEFS = -D__PGI -D__FFTW

9
src/MACHINEDEFS.SGI
View file

@ -6,15 +6,16 @@ CPP = /usr/lib/cpp
FC = f90 -freeform
FCfixed = f90 -fixedform
LDR = f90
AR = ar -r
CDEFS = -D__IRIX -D__FFTW
CDEFS = -D__IRIX -D__FFTSG -D__FFTW
DEBUG = -g -C
OPT = -O3 -Ofast
BFLAGS = -u $(CDEFS) -automatic -64
OPT = -O3
BFLAGS = -u $(CDEFS) -automatic
FFLAGS = $(BFLAGS) $(OPT) -macro_expand
LFLAGS = $(BFLAGS) $(OPT)
LIBS = -L${HOME}/lib -lfftw -lcomplib.sgimath
LIBS = -L/usr/local/lib -lfftw -lcomplib.sgimath
OBJECTS_ARCHITECTURE = machine_irix.o

1
src/MACHINEDEFS.T3E
View file

@ -6,6 +6,7 @@ CPP = cpp
FC = f90 -f free
FCfixed = f90 -f fixed
LDR = f90
AR = ar -r
CDEFS = -D__T3E -D__FFTW

11
src/Makefile
View file

@ -1,7 +1,8 @@
.SUFFIXES: .o .F .d
CP2KHOME= ${HOME}/CP2K
CP2KHOME= ${HOME}/CP2K/
FORPAR = $(CP2KHOME)/tools/forpar.x -chkint
SFMAKEDEPEND = $(CP2KHOME)/tools/sfmakedepend -m int -s -f
include OBJECTDEFS
@ -9,6 +10,8 @@ include MACHINEDEFS
OBJECTS = $(OBJECTS_GENERIC) $(OBJECTS_ARCHITECTURE)
LIBRARIES = $(LIBS) -L./lib -lfftsg
#################################
PROG = cp2k.x
@ -16,7 +19,7 @@ PROG = cp2k.x
all: $(PROG)
$(PROG): $(OBJECTS)
$(LDR) $(LFLAGS) -o $(PROG) $(OBJECTS) $(LIBS)
$(LDR) $(LFLAGS) -o $(PROG) $(OBJECTS) $(LIBRARIES)
%.o : %.F
$(FC) -c $(FFLAGS) $*.F
@ -26,13 +29,13 @@ $(PROG): $(OBJECTS)
parallel_include.o : parallel_include.F
$(FCfixed) -c $(FFLAGS) $*.F
$(CPP) -P $(CDEFS) -U__parallel $*.F $*.f
$(CPP) -P $(CDEFS) $*.F $*.f
-$(FORPAR) -fix $*.f
rm -f $*.f
%.d : %.F
$(CPP) -P $(CDEFS) $*.F $*.f
$(PERL) $(CP2KHOME)/tools/sfmakedepend -m int -s -f $*.d $*.f
$(PERL) $(SFMAKEDEPEND) $*.d $*.f
rm -f $*.f $*.d.old
depend:

16
src/OBJECTDEFS
View file

@ -1,16 +1,20 @@
OBJECTS_GENERIC = \
atoms_input.o coefficients.o coefficient_types.o \
constraint.o convert_units.o dgs.o dg_rho0s.o dg_types.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 \
fft_tools.o fist.o fist_debug.o fist_global.o fist_input.o \
fermi.o fft_tools.o fftsg_lib.o fftw_lib.o \
fist_debug.o fist_global.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 \
kinds.o linklists.o linklist_control.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_utilities.o linklist_verlet_list.o machine.o \
mathconstants.o mol_force.o molecule_types.o \
@ -21,5 +25,5 @@ physcon.o pme.o pws.o pw_types.o pw_grids.o \
pw_grid_types.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.o tbmd_debug.o tbmd_force.o tbmd_global.o tbmd_input.o \
tbmd_types.o timings.o unit.o util.o
tbmd_debug.o tbmd_force.o tbmd_global.o tbmd_initialize.o \
tbmd_input.o tbmd_types.o timings.o unit.o util.o

318
src/band.F Normal file
View file

@ -0,0 +1,318 @@
!------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!------------------------------------------------------------------------------!
!
MODULE band
!
!------------------------------------------------------------------------------!
USE kinds, ONLY : dbl
USE stop_program, ONLY : stop_prg, stop_memory
USE brillouin, ONLY : kpoint_type
USE fermi, ONLY : fermi_distribution_type
!
IMPLICIT NONE
!
PRIVATE
!
PUBLIC :: band_structure_type, init_band_structure, band_structure_info, &
adjust_mu, occupation
!
TYPE band_structure_type
TYPE (fermi_distribution_type), POINTER :: fd
TYPE (kpoint_type), POINTER :: kpt
REAL (dbl), DIMENSION (:,:,:), POINTER :: ev
REAL (dbl), DIMENSION (:,:,:), POINTER :: oc
END TYPE band_structure_type
!
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
!
!------------------------------------------------------------------------------!
SUBROUTINE init_band_structure(bs,fd,kp)
IMPLICIT NONE
TYPE (band_structure_type), INTENT (INOUT) :: bs
TYPE (fermi_distribution_type), INTENT (IN), TARGET :: fd
TYPE (kpoint_type), INTENT (IN), TARGET :: kp
INTEGER :: nspin, ns, nk, isos
bs%fd => fd
bs%kpt => kp
nspin = fd%spin_polarization + 1
ns = max(fd%na,fd%nb)
nk = kp%nkpt
IF (associated(bs%ev)) NULLIFY (bs%ev)
IF (associated(bs%oc)) NULLIFY (bs%oc)
ALLOCATE (bs%ev(ns,nk,nspin),STAT=isos)
IF (isos/=0) CALL stop_memory('init_band_structure', &
'bs%ev',ns*nk*nspin)
ALLOCATE (bs%oc(ns,nk,nspin),STAT=isos)
IF (isos/=0) CALL stop_memory('init_band_structure', &
'bs%oc',ns*nk*nspin)
END SUBROUTINE init_band_structure
!------------------------------------------------------------------------------!
SUBROUTINE band_structure_info(bs,punit)
IMPLICIT NONE
TYPE (band_structure_type), INTENT (IN) :: bs
INTEGER, INTENT (IN) :: punit
INTEGER :: i, ik, left, nm, nk
LOGICAL :: nokp
WRITE (punit,'(/,1x,79("-"))')
WRITE (punit,'(" -",77x,"-")')
WRITE (punit,'(" -",25x,a,24x,"-")') ' B A N D S T R U C T U R E '
WRITE (punit,'(" -",77x,"-")')
WRITE (punit,'(1x,79("-"))')
WRITE (punit,'(" -",T23,a,f8.3,T80,"-")') &
' Chemical Potential [a.u.] = ', bs%fd%mu
WRITE (punit,'(1x,79("-"))')
nokp = bs%kpt%scheme=='GAMMA' .OR. bs%kpt%scheme=='NULL'
IF (nokp) THEN
nk = 1
ELSE
nk = bs%kpt%nkpt
END IF
DO ik = 1, nk
IF (nokp) THEN
WRITE (punit,'(1x,79("."))')
WRITE (punit, &
'(" .",a,i5,6x,a,f8.3,a,f8.3,a,f8.3,4x,a,f8.5,T80,".")') &
' K-point:', ik, ' x=', bs%kpt%xk(1,ik), ' y=', bs%kpt%xk(2,ik), &
' z=', bs%kpt%xk(3,ik), ' weight=', bs%kpt%weight(ik)
WRITE (punit,'(1x,79("."))')
END IF
IF (bs%fd%spin_polarization==0) THEN
WRITE (punit,'(A,A)') ' State Occupation Eigenvalue[au] ', &
' State Occupation Eigenvalue[au] '
DO i = 1, bs%fd%nstate, 2
left = bs%fd%nstate - i + 1
IF (left>1) THEN
WRITE (punit,'(I5,5x,F10.3,6x,F13.6,T42,I5,5x,F10.3,6x,F13.6)' &
) i, bs%oc(i,ik,1), bs%ev(i,ik,1), i + 1, bs%oc(i+1,ik,1), &
bs%ev(i+1,ik,1)
ELSE
WRITE (punit,'(I5,5x,F10.3,6x,F13.6)') i, bs%oc(i,ik,1), &
bs%ev(i,ik,1)
END IF
END DO
ELSE
WRITE (punit,'(T12,A,T54,A)') ' Alpha Electrons ', &
' Beta Electrons'
WRITE (punit,'(A,A)') ' State Occupation Eigenvalue[au] ', &
' State Occupation Eigenvalue[au] '
nm = min(bs%fd%na,bs%fd%nb)
DO i = 1, nm
WRITE (punit,'(I5,5x,F10.3,6x,F13.6,T42,I5,5x,F10.3,6x,F13.6)') &
i, bs%oc(i,ik,1), bs%ev(i,ik,1), i, bs%oc(i,ik,2), &
bs%ev(i,ik,2)
END DO
DO i = nm + 1, bs%fd%na
WRITE (punit,'(I5,5x,F10.3,6x,F13.6)') i, bs%oc(i,ik,1), &
bs%ev(i,ik,1)
END DO
DO i = nm + 1, bs%fd%nb
WRITE (punit,'(T40,I5,5x,F10.3,6x,F13.6)') i, bs%oc(i,ik,2), &
bs%ev(i,ik,2)
END DO
END IF
END DO
WRITE (punit,'(1x,79("-"),/)')
END SUBROUTINE band_structure_info
!------------------------------------------------------------------------------!
SUBROUTINE adjust_mu(bs)
!
! adjust the chemical potential with the bisection method
! eigenvalues do not have to be sorted
! routine adapted from CPMD
!
IMPLICIT NONE
TYPE (band_structure_type), INTENT (INOUT) :: bs
REAL (dbl) :: amu, amu1, amu2, damu, sdeg
REAL (dbl) :: rhint, rhint1, rhint2, precision
INTEGER, PARAMETER :: it_max = 200
INTEGER :: nk, it
!
precision = precpar*log(float(bs%fd%nel+1))
!..min and max eigenvalue
amu1 = minval(bs%ev)
amu2 = maxval(bs%ev)
!
DO
rhint1 = rhoint(bs,amu1)
IF (rhint1<0._dbl) EXIT
amu1 = -2._dbl*abs(amu1)
END DO
rhint2 = rhoint(bs,amu2)
IF (rhint2<0._dbl) CALL stop_prg('ADJUST_MU', &
'failed to find chem. potential','number of states is too small')
damu = amu2 - amu1
amu = amu1
!
DO it = 1, it_max
amu = amu1 + 0.5_dbl*damu
rhint = rhoint(bs,amu)
IF (damu<delta) THEN
IF (abs(rhint)<precision) THEN
EXIT
ELSE IF (abs(rhint1)<precision) THEN
amu = amu1
EXIT
ELSE IF (abs(rhint1)<precision) THEN
amu = amu2
EXIT
END IF
END IF
IF (rhint<0._dbl) THEN
amu1 = amu
rhint1 = rhint
ELSE
amu2 = amu
rhint2 = rhint
END IF
damu = amu2 - amu1
END DO
IF (damu>delta) THEN
CALL stop_prg('ADJUST_MU','failed to find chem. potential')
END IF
bs%fd%mu = amu
!
END SUBROUTINE adjust_mu
!------------------------------------------------------------------------------!
FUNCTION rhoint(bs,mu)
IMPLICIT NONE
TYPE (band_structure_type), INTENT (IN) :: bs
REAL (dbl), INTENT (IN) :: mu
REAL (dbl) :: rhoint, xlm, arg, argmin, argmax, sdeg, remain
REAL (dbl) :: dramu, drwe, wbigwe
INTEGER :: is, ik, i, n
argmin = log(toll)
argmax = -argmin
IF (bs%fd%spin_polarization==0) THEN
sdeg = 2._dbl
ELSE
sdeg = 1._dbl
END IF
dramu = dround(mu,delta)
rhoint = 0._dbl
wbigwe = 0._dbl
DO is = 1, size(bs%ev(1,1,:))
DO ik = 1, size(bs%ev(1,:,1))
IF (is==1) n = bs%fd%na
IF (is==2) n = bs%fd%nb
DO i = 1, n
arg = -bs%fd%betael*(bs%ev(i,ik,is)-mu)
IF (arg>argmax) THEN
xlm = 1._dbl
ELSE IF (arg<argmin) THEN
xlm = 0._dbl
ELSE
xlm = exp(arg)/(exp(arg)+1._dbl)
END IF
drwe = dround(bs%ev(i,ik,is),delta)
IF (abs(drwe-dramu)<delta) THEN
wbigwe = wbigwe + bs%kpt%weight(ik)
ELSE
rhoint = rhoint + bs%kpt%weight(ik)*xlm
END IF
END DO
END DO
END DO
remain = float(bs%fd%nel) - rhoint
IF (remain<0._dbl) remain = sdeg
IF (wbigwe>0._dbl) THEN
remain = remain/wbigwe
remain = min(remain,sdeg)
DO is = 1, size(bs%ev(1,1,:))
DO ik = 1, size(bs%ev(1,:,1))
IF (is==1) n = bs%fd%na
IF (is==2) n = bs%fd%nb
DO i = 1, n
drwe = dround(bs%ev(i,ik,is),delta)
IF (abs(drwe-dramu)<delta) rhoint = rhoint + &
bs%kpt%weight(ik)*remain
END DO
END DO
END DO
END IF
rhoint = sdeg*rhoint - float(bs%fd%nel)
END FUNCTION rhoint
!------------------------------------------------------------------------------!
SUBROUTINE occupation(bs)
IMPLICIT NONE
TYPE (band_structure_type), INTENT (INOUT) :: bs
REAL (dbl) :: xlm, arg, argmin, argmax, sdeg
REAL (dbl) :: dramu, drwe, wbigwe, remain, rhotot
INTEGER :: is, ik, i, n
argmin = log(toll)
argmax = -argmin
IF (bs%fd%spin_polarization==0) THEN
sdeg = 2._dbl
ELSE
sdeg = 1._dbl
END IF
dramu = dround(bs%fd%mu,delta)
wbigwe = 0._dbl
rhotot = 0._dbl
DO is = 1, size(bs%oc(1,1,:))
DO ik = 1, size(bs%oc(1,:,1))
IF (is==1) n = bs%fd%na
IF (is==2) n = bs%fd%nb
DO i = 1, n
arg = -bs%fd%betael*(bs%ev(i,ik,is)-bs%fd%mu)
IF (arg>argmax) THEN
xlm = 1._dbl
ELSE IF (arg<argmin) THEN
xlm = 0._dbl
ELSE
xlm = exp(arg)/(exp(arg)+1._dbl)
END IF
drwe = dround(bs%ev(i,ik,is),delta)
IF (abs(drwe-dramu)<delta) THEN
wbigwe = wbigwe + bs%kpt%weight(ik)
ELSE
rhotot = rhotot + bs%kpt%weight(ik)*xlm
bs%oc(i,ik,is) = xlm*sdeg
END IF
END DO
END DO
END DO
remain = float(bs%fd%nel) - rhotot
IF (remain<0._dbl) remain = sdeg
IF (wbigwe>0._dbl) THEN
remain = remain/wbigwe
remain = min(remain,sdeg)
DO is = 1, size(bs%ev(1,1,:))
DO ik = 1, size(bs%ev(1,:,1))
IF (is==1) n = bs%fd%na
IF (is==2) n = bs%fd%nb
DO i = 1, n
drwe = dround(bs%ev(i,ik,is),delta)
IF (abs(drwe-dramu)<delta) bs%oc(i,ik,is) = remain
END DO
END DO
END DO
END IF
END SUBROUTINE occupation
!------------------------------------------------------------------------------!
FUNCTION dround(a,delta)
IMPLICIT NONE
REAL (dbl), INTENT (IN) :: a, delta
REAL (dbl) :: dround
dround = a - mod(a,delta)
END FUNCTION dround
!------------------------------------------------------------------------------!
END MODULE band
!------------------------------------------------------------------------------!

203
src/brillouin.F Normal file
View file

@ -0,0 +1,203 @@
!------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!------------------------------------------------------------------------------!
!
MODULE brillouin
!
!------------------------------------------------------------------------------!
USE kinds, ONLY : dbl
USE global_types, ONLY : global_environment_type
USE stop_program, ONLY : stop_prg, stop_memory
USE string_utilities, ONLY : uppercase
USE parser, ONLY : parser_init, parser_end, read_line, test_next, &
cfield, p_error, get_real, get_int, stop_parser
!
IMPLICIT NONE
!
PRIVATE
!
PUBLIC :: kpoint_type, brillouin_info, kpoint_input
!
TYPE kpoint_type
CHARACTER (len=20) :: scheme
INTEGER :: nk(3)
REAL (dbl) :: shift(3)
LOGICAL :: symmetry
INTEGER :: wfn_type
INTEGER :: nkpt
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] !
!! !
!!<----------------------------------------------------------------------------!
SUBROUTINE kpoint_input(kp,inpar)
IMPLICIT NONE
TYPE (kpoint_type), INTENT (OUT) :: kp
TYPE (global_environment_type), INTENT (IN) :: inpar
CHARACTER (len=20) :: string, str2
CHARACTER (len=6) :: label
INTEGER :: iw, ierror, ilen, isos, i, source, group
!..defaults
kp%scheme = 'NULL'
kp%symmetry = .FALSE.
kp%wfn_type = 0
kp%nkpt = 0
kp%nk = 0
kp%shift = 0._dbl
iw = inpar%scr
!..parse the input section
label = '&KPOINT'
CALL parser_init(inpar%input_file_name,label,ierror,inpar)
IF (ierror/=0) THEN
IF (inpar%ionode) &
WRITE (iw,'(a,a)') ' No input section &KPOINT found on file ', &
adjustl(inpar%input_file_name)
ELSE
CALL read_line
DO WHILE (test_next()/='X')
ilen = 8
CALL cfield(string,ilen)
CALL uppercase(string)
SELECT CASE (string)
CASE DEFAULT
CALL p_error()
CALL stop_parser('kpoint_input','unknown option')
CASE ('SCHEME')
ilen = 20
CALL cfield(str2,ilen)
CALL uppercase(str2)
SELECT CASE (str2)
CASE DEFAULT
CALL p_error()
CALL stop_prg('kpoint_input','Scheme: unknown option')
CASE ('GAMMA')
kp%scheme = str2
CASE ('MONKHORST-PACK')
kp%scheme = str2
kp%nk(1) = get_int()
kp%nk(2) = get_int()
kp%nk(3) = get_int()
CASE ('MACDONALD')
kp%scheme = str2
kp%nk(1) = get_int()
kp%nk(2) = get_int()
kp%nk(3) = get_int()
kp%shift(1) = get_real()
kp%shift(2) = get_real()
kp%shift(3) = get_real()
CASE ('GENERAL')
kp%scheme = str2
kp%nkpt = get_int()
ALLOCATE (kp%xk(3,kp%nkpt),STAT=isos)
IF (isos/=0) CALL stop_memory('kpoint_input', &
'kp%xk',3*kp%nkpt)
ALLOCATE (kp%weight(kp%nkpt),STAT=isos)
IF (isos/=0) CALL stop_memory('kpoint_input', &
'kp%weight',kp%nkpt)
DO i = 1, kp%nkpt
kp%xk(1,i) = get_real()
kp%xk(2,i) = get_real()
kp%xk(3,i) = get_real()
kp%weight(i) = get_real()
END DO
END SELECT
CASE ('SYMMETRY')
ilen = 3
CALL cfield(str2,ilen)
CALL uppercase(str2)
IF (str2=='OFF') kp%symmetry = .FALSE.
IF (str2=='ON') kp%symmetry = .TRUE.
CASE ('WAVEFUNC')
ilen = 4
CALL cfield(str2,ilen)
CALL uppercase(str2)
SELECT CASE (str2)
CASE DEFAULT
CALL p_error()
CALL stop_parser('kpoint_input',&
'Wavefunctions: unknown option')
CASE ('REAL')
kp%wfn_type = 0
CASE ('COMP')
kp%wfn_type = 1
END SELECT
END SELECT
END DO
END IF
CALL parser_end
!..update defaults
IF (kp%scheme=='NULL' .OR. kp%scheme=='GAMMA') THEN
kp%nkpt = 1
ALLOCATE (kp%xk(3,1),STAT=isos)
IF (isos/=0) CALL stop_memory('kpoint_input','kp%xk',3)
ALLOCATE (kp%weight(1),STAT=isos)
IF (isos/=0) CALL stop_memory('kpoint_input','kp%weight',1)
kp%xk(:,:) = 0._dbl
kp%weight(1) = 1._dbl
kp%scheme = 'GAMMA'
END IF
END SUBROUTINE kpoint_input
!------------------------------------------------------------------------------!
SUBROUTINE brillouin_info(kp,punit)
IMPLICIT NONE
TYPE (kpoint_type), INTENT (IN) :: kp
INTEGER, INTENT (IN) :: punit
INTEGER :: i
IF (kp%scheme=='GAMMA') THEN
WRITE (punit,*)
WRITE (punit,'(A,T57,A)') ' BRILLOUIN|', ' Gamma-point calculation'
WRITE (punit,'(A,T76,A)') ' BRILLOUIN| Wavefunction type', ' REAL'
ELSE
WRITE (punit,*)
WRITE (punit,'(A,T61,A)') ' BRILLOUIN| K-point scheme ', &
adjustr(kp%scheme)
IF (kp%scheme=='MONKHORST-PACK') THEN
WRITE (punit,'(A,T66,3I5)') ' BRILLOUIN| K-Point grid', kp%nk
ELSE IF (kp%scheme=='MACDONALD') THEN
WRITE (punit,'(A,T66,3I5)') ' BRILLOUIN| K-Point grid', kp%nk
WRITE (punit,'(A,T51,3F10.4)') ' BRILLOUIN| K-Point shift', &
kp%shift
END IF
IF (kp%symmetry) THEN
WRITE (punit,'(A,T76,A)') ' BRILLOUIN| K-Point symmetry', ' ON'
ELSE
WRITE (punit,'(A,T76,A)') ' BRILLOUIN| K-Point symmetry', ' OFF'
END IF
IF (kp%wfn_type==0) THEN
WRITE (punit,'(A,T76,A)') ' BRILLOUIN| Wavefunction type', ' REAL'
ELSE
WRITE (punit,'(A,T73,A)') ' BRILLOUIN| Wavefunction type', &
' COMPLEX'
END IF
WRITE (punit,'(A,T71,I10)') ' BRILLOUIN| Number of K-points ', &
kp%nkpt
WRITE (punit,'(A,T30,A,T48,A,T63,A,T78,A)') ' BRILLOUIN| Number ', &
'Weight', 'X', 'Y', 'Z'
DO i = 1, kp%nkpt
WRITE (punit,'(A,I5,3X,4F15.5)') ' BRILLOUIN| ', i, kp%weight(i), &
kp%xk(1,i), kp%xk(2,i), kp%xk(3,i)
END DO
END IF
END SUBROUTINE brillouin_info
!------------------------------------------------------------------------------!
END MODULE brillouin
!------------------------------------------------------------------------------!

178
src/cntl_input.F Normal file
View file

@ -0,0 +1,178 @@
!-----------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
MODULE cntl_input
USE ewald_parameters_types, ONLY : ewald_parameters_type
USE global_types, ONLY : global_environment_type
USE input_types, ONLY : setup_parameters_type
USE kinds, ONLY : dbl
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, xstring
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]
!! !
!!<---------------------------------------------------------------------------!
SUBROUTINE read_cntl_section ( setup, ewald_param, globenv )
IMPLICIT NONE
! Arguments
TYPE ( setup_parameters_type ), INTENT ( INOUT ) :: setup
TYPE ( ewald_parameters_type ), INTENT ( INOUT ) :: ewald_param
TYPE ( global_environment_type ), INTENT ( INOUT ) :: globenv
! Locals
INTEGER :: ierror, ilen, ia, ie, i, j, n, iw, source, group
CHARACTER ( LEN = 20 ) :: string, str2
CHARACTER ( LEN = 5 ) :: label
CHARACTER ( LEN = 3 ), PARAMETER :: yn ( 0:1 ) = (/ ' NO', 'YES' /)
!------------------------------------------------------------------------------
!..defaults
setup % run_type = 'MD'
setup % unit_type = 'KELVIN'
setup % perd = 1
setup % symmetry = .false.
ewald_param % alpha = 0.4_dbl
ewald_param % gmax = 10
ewald_param % ns_max = 10
ewald_param % epsilon = 1.e-6_dbl
ewald_param % ewald_type = 'NONE'
CALL xstring(globenv % project_name,ia,ie)
setup % set_file_name = globenv % project_name(ia:ie) // '.set'
setup % input_file_name = globenv % project_name(ia:ie) // '.dat'
iw = globenv % scr
!..parse the input section
label = '&CNTL'
CALL parser_init(globenv % input_file_name,label,ierror,globenv)
IF (ierror /= 0 ) THEN
IF (globenv % ionode) &
WRITE ( iw, '( a )' ) ' No input section &CNTL found '
ELSE
CALL read_line
DO WHILE (test_next()/='X')
ilen = 8
CALL cfield ( string, ilen )
CALL uppercase ( string )
SELECT CASE ( string )
CASE DEFAULT
CALL p_error()
CALL stop_parser( 'read_cntl_section','unknown option')
CASE ( 'SIMULATI')
ilen = 20
CALL cfield(setup % run_type,ilen)
CALL uppercase(setup % run_type )
CASE ( 'PRINTLEV')
globenv % print_level = get_int()
CASE ( 'UNITS')
ilen = 20
CALL cfield(setup % unit_type,ilen)
CALL uppercase(setup % unit_type )
CASE ( 'PERIODIC')
setup % perd(1) = get_int()
setup % perd(2) = get_int()
setup % perd(3) = get_int()
CASE ( 'EWALD_TY')
ilen=20
CALL cfield(string,ILEN)
CALL uppercase ( string )
SELECT CASE(string)
CASE( 'EWALD_GAUSS')
ewald_param % ewald_type = 'ewald_gauss'
CALL uppercase(ewald_param % ewald_type )
CASE( 'PME_GAUSS')
ewald_param % ewald_type = 'pme_gauss'
CALL uppercase(ewald_param % ewald_type )
END SELECT
! if no type specified, assume ewald_gauss
CASE ( 'EWALD_PA')
ewald_param % alpha = get_real()
SELECT CASE (ewald_param % ewald_TYPE ( 1:3))
CASE DEFAULT
ewald_param % gmax = get_int()
CASE ( 'PME')
ewald_param % ns_max = get_int()
IF ( test_next() == 'N' ) THEN
ewald_param % epsilon = get_real()
END IF
END SELECT
CASE ('SYMMETRY')
ilen = 3
CALL cfield(str2,ilen)
CALL uppercase(str2)
IF ( str2(1:2) == "ON" ) setup%symmetry=.true.
IF ( str2(1:3) == "OFF" ) setup%symmetry=.false.
CASE ( 'SET_FILE')
ilen = 20
CALL cfield(setup % set_file_name,ilen)
CASE ( 'INPUT_FI')
ilen = 20
CALL cfield(setup % input_file_name,ilen)
END SELECT
! check for trailing rubbish
CALL read_line
END DO
END IF
CALL parser_end
!..end of parsing the input section
!..write some information to output
IF (globenv % ionode) THEN
IF ( globenv % print_level >= 0 ) THEN
WRITE ( iw, '( A, T71, A )' ) &
' CONTROL| Run type ', ADJUSTR ( setup % run_type )
WRITE ( iw, '( A, T71, A )' ) &
' CONTROL| Unit type ', ADJUSTR ( setup % unit_type )
WRITE ( iw, '( A, T78, A )' ) ' CONTROL| Periodic in X direction ', &
yn(setup % perd(1))
WRITE ( iw, '( A, T78, A )' ) ' CONTROL| Periodic in Y direction ', &
yn(setup % perd(2))
WRITE ( iw, '( A, T78, A )' ) ' CONTROL| Periodic in Z direction ', &
yn(setup % perd(3))
IF ( setup%symmetry ) THEN
WRITE (iw,'(A,T78,A)') ' CONTROL| Use of symmetry ','Yes'
ELSE
WRITE (iw,'(A,T78,A)') ' CONTROL| Use of symmetry ',' No'
END IF
WRITE ( iw, '( A, T61, A )' ) ' CONTROL| Set file name', &
ADJUSTR ( setup % set_file_name )
WRITE ( iw, '( A, T61, A )' ) ' CONTROL| Input file name', &
ADJUSTR ( setup % input_file_name )
WRITE ( iw, '( A, T76, I5 )' ) &
' CONTROL| Print level ', globenv % print_level
WRITE ( iw, '( )' )
END IF
END IF
END SUBROUTINE read_cntl_section
!******************************************************************************
END MODULE cntl_input

413
src/control_module.F Normal file
View file

@ -0,0 +1,413 @@
!-----------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
MODULE control_module
USE atoms_input, ONLY : read_coord_vel, system_type
USE convert_units, ONLY : convert
USE global_types, ONLY : global_environment_type
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, tbmd_header
USE kinds, ONLY : dbl
USE band, ONLY : band_structure_type, init_band_structure
USE brillouin, ONLY : kpoint_type, brillouin_info, kpoint_input
USE kpoint_initialization, ONLY : initialize_kpoints
USE fermi, ONLY : fermi_distribution_type, init_fermi_dist, fermi_info
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
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, &
charge
USE molecule_types, ONLY : molecule_type, intra_parameters_type
USE nose, ONLY : extended_parameters_type
USE pair_potential, ONLY : 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, interaction_type
USE timings, ONLY : timeset, timestop, trace_debug
USE unit, ONLY : unit_convert_type, set_units
USE util, ONLY : close_unit, get_share
IMPLICIT NONE
PRIVATE
PUBLIC :: control
TYPE ( mdio_parameters_type ) :: mdio
CONTAINS
!-----------------------------------------------------------------------------!
! CNTL CNTL CNTL CNTL CNTL CNTL CNTL CNTL CNTL CNTL CNTL CNTL CNTL CNTL CNTL !
!-----------------------------------------------------------------------------!
SUBROUTINE control ( globenv )
IMPLICIT NONE
TYPE ( global_environment_type ), INTENT ( INOUT ) :: globenv
! Locals
REAL ( dbl ) :: cons, ecut, qi, qj
REAL ( dbl ), DIMENSION ( :, : ), ALLOCATABLE :: rcut
INTEGER :: handle1, handle2, handle3 , itimes, isos, para_comm_cart
INTEGER :: iat, jat, nel
CHARACTER ( LEN = ATOMNAMESLENGTH ), DIMENSION ( : ), POINTER :: atom_names
CHARACTER ( LEN = 20 ) :: set_fn
LOGICAL :: tbmd, fist
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 ( interaction_type ) :: inter
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
TYPE ( kpoint_type ) :: kp
TYPE ( fermi_distribution_type ) :: fd
TYPE ( band_structure_type ) :: bs
!------------------------------------------------------------------------------
! IF( globenv % ionode ) CALL trace_debug ( "START" )
CALL timeset ( 'CONTROL', 'I', ' ', handle1 )
CALL timeset ( 'CNTL_INIT', 'I', ' ', handle2 )
IF ( globenv % program_name == 'FIST' ) THEN
fist = .true.
tbmd = .false.
ELSE IF ( globenv % program_name == 'TBMD' ) THEN
tbmd = .true.
fist = .false.
ELSE
call stop_prg ( ' control ',' program_name not specified' )
ENDIF
IF ( globenv % ionode ) THEN
IF ( fist ) CALL fist_header ( globenv % scr )
IF ( tbmd ) CALL tbmd_header ( globenv % scr )
END IF
!..read cntl section
CALL read_cntl_section ( setup, ewald_param, globenv )
! read from the setup and molecule section of *.set
set_fn = setup % set_file_name
CALL read_setup_section ( mol_setup, set_fn, globenv )
CALL read_molecule_section ( mol_setup, set_fn, globenv )
! read force_field information for classical MD
! read pair potential information for TB
CALL read_force_field_section ( setup, mol_setup, set_fn, &
intra_param, inter%potparm, atom_names, pstat, globenv )
!..read Hamiltonian section
IF ( tbmd ) THEN
CALL read_tb_hamiltonian ( setup, atom_names, inter%tbatom, globenv )
CALL read_tb_hopping_elements ( setup, atom_names, &
inter%tbatom, inter%tbhop, globenv )
END IF
!..read the input of the molecular dynamics section
CALL read_md_section ( simpar, globenv, mdio )
simpar % program = globenv % program_name
!..read atomic coordinates, velocities (optional) and the simulation box
ainp % rtype = simpar % read_type
CALL read_coord_vel ( ainp, setup % input_file_name, globenv )
!..initialize box, perd
struc % box % hmat = ainp % box
struc % box % perd = setup % perd
!..K-Points
IF ( SUM ( struc % box % perd ) == 0 ) THEN
kp%scheme = "NULL"
ELSE IF ( tbmd ) THEN
CALL kpoint_input ( kp, globenv )
ELSE
kp%scheme = "NULL"
END IF
!..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, globenv )
!..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 = inter%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, globenv )
!..initialize extended_parameters_type and get number of degrees of freedom
CALL initialize_extended_type ( struc % box, simpar, &
struc % molecule, mol_setup, nhcp, globenv )
! initialize velocities if needed
IF ( simpar % read_type == 'POS' ) THEN
CALL initialize_velocities ( simpar, struc % part, globenv )
END IF
!...initialize splines
inter%potparm ( :, : ) % energy_cutoff = 0.0_dbl
inter%potparm ( :, : ) % e_cutoff_coul = 0.0_dbl
CALL spline_nonbond_control ( inter%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 ( 'control', 'rcut', 0 )
rcut ( :, : ) = inter%potparm ( :, : ) % rcutsq
CALL set_ll_parm ( globenv, simpar % verlet_skin, &
setup % natom_type, rcut, simpar % n_cell )
CALL set_ll_parm ( globenv, printlevel = globenv % print_level, &
ltype = 'NONBOND' )
DEALLOCATE ( rcut, STAT = isos )
IF ( isos /= 0 ) CALL stop_memory ( 'control', 'rcut' )
!..deallocate arrays needed for atom input
IF ( ASSOCIATED ( ainp % c ) ) THEN
DEALLOCATE ( ainp % c, STAT = isos )
IF ( isos /= 0 ) CALL stop_memory ( 'control', 'ainp%c' )
END IF
IF ( ASSOCIATED ( ainp % v ) ) THEN
DEALLOCATE ( ainp % v, STAT = isos )
IF ( isos /= 0 ) CALL stop_memory ( 'control', 'ainp%v' )
END IF
!
!..initialize the on-site terms for TB
IF ( tbmd ) CALL tbmd_init ( struc % part, tbatom, tbhop )
!
!..Symmetry and K-points
!
IF ( SUM ( struc % box % perd ) == 0 ) THEN
!!!!!! symmetry setup for molecules
ELSE IF ( tbmd ) THEN
CALL initialize_kpoints ( globenv, kp, setup % symmetry, &
struc % box % hmat, struc % part )
IF ( globenv%ionode .AND. globenv%print_level > 0) &
CALL brillouin_info ( kp, globenv%scr )
END IF
!..nota bene: we use the defaults for the elctron temperature and
! spin polarisation
IF ( tbmd ) then
nel = tb_get_numel ( 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)
ENDIF
CALL timestop ( zero, handle2 )
CALL timeset ( 'CNTL_WORK', 'I', ' ', handle3 )
IF ( setup % run_type == 'DEBUG' ) THEN
!..debug the forces
!..initialize integrator
CALL set_integrator ( globenv, mdio )
IF ( fist ) then
CALL fist_debug_control( globenv, ewald_param, struc%part, struc%pnode, &
struc%molecule, struc%box, thermo, inter%potparm, simpar % ensemble )
ELSE IF ( tbmd ) THEN
CALL tbmd_debug_control( globenv, ewald_param, struc%part, struc%pnode, &
struc%molecule, struc%box, thermo, inter%potparm, simpar % ensemble )
ENDIF
ELSE
!..initialize integrator
CALL set_integrator ( globenv, mdio )
!..MD
itimes = 0
CALL force ( struc, inter, thermo, simpar, ewald_param, globenv )
IF ( globenv % ionode .AND. ewald_param % ewald_type /= 'NONE' ) &
CALL ewald_print ( globenv % 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, inter, thermo, &
struc, ewald_param, nhcp )
IF ( MOD ( itimes, mdio % idump ) == 0 ) &
CALL dump_variables ( struc, mdio % dump_file_name, globenv )
END DO
CALL dump_variables ( struc, mdio % dump_file_name, globenv )
IF ( globenv % ionode ) CALL close_unit ( 10, 99 )
END IF
!..deallocate memory for atoms and molecules
CALL deallocmem ( struc )
CALL timestop ( zero, handle3 )
CALL timestop ( zero, handle1 )
END SUBROUTINE control
!-----------------------------------------------------------------------------!
! CNTL CNTL CNTL CNTL CNTL CNTL CNTL CNTL CNTL CNTL CNTL CNTL CNTL CNTL CNTL !
!-----------------------------------------------------------------------------!
SUBROUTINE allocmem ( ainp, mol_setup, struc, globenv )
IMPLICIT NONE
! Arguments
TYPE ( molecule_type ), DIMENSION ( : ), INTENT ( IN ) :: mol_setup
TYPE ( structure_type ), INTENT ( INOUT ) :: struc
TYPE ( system_type ), INTENT ( IN ) :: ainp
TYPE ( global_environment_type ), INTENT ( IN ) :: globenv
! Locals
INTEGER :: iw, natoms, nnodes, nmol, nmoltype, ios, iat, i, nsh
!------------------------------------------------------------------------------
struc % name = globenv % program_name // ' MOLECULAR SYSTEM'
IF ( globenv % num_pe == 1 ) THEN
natoms = SIZE ( ainp % c ( 1, : ) )
ALLOCATE ( struc % part ( 1:natoms ), STAT = ios )
IF ( ios /= 0 ) CALL stop_memory ( 'control', 'part', natoms )
ALLOCATE ( struc % pnode ( 1:natoms ), STAT = ios )
IF ( ios /= 0 ) CALL stop_memory ( 'control', 'pnode', natoms )
nmol = SUM ( mol_setup ( : ) % num_mol )
ALLOCATE ( struc % molecule ( 1:nmol ), STAT = ios )
IF ( ios /= 0 ) CALL stop_memory ( 'control', 'molecule', nmol )
IF ( globenv % ionode .AND. globenv % print_level > 3 ) THEN
iw = globenv % scr
WRITE ( iw, '( A )' )
WRITE ( iw, '( A, T71, I10 )' ) &
' CONTROL| Number of allocated particles ', natoms
WRITE ( iw,'( A, T71, I10 )' ) &
' CONTROL| Number of allocated particle nodes ', natoms
WRITE ( iw, '( A, T71, I10 )' ) &
' CONTROL| 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 ( 'control', 'part', natoms )
nmoltype = SIZE ( mol_setup )
nmol = 0
nnodes = 0
DO i = 1, nmoltype
nsh = get_share ( mol_setup ( i ) % num_mol, &
globenv % num_pe, globenv % mepos )
nmol = nmol + nsh
nnodes = nnodes + nsh * mol_setup ( i ) % molpar % natom
END DO
ALLOCATE ( struc % molecule ( 1:nmol ), STAT = ios )
IF ( ios /= 0 ) CALL stop_memory ( 'control', 'molecule' , nmol )
ALLOCATE ( struc % pnode ( 1:nnodes ), STAT = ios )
IF ( ios /= 0 ) CALL stop_memory ( 'control', 'pnode', nnodes )
IF ( globenv % ionode .AND. globenv % print_level > 3 ) THEN
iw = globenv % scr
WRITE ( iw, '( A )' )
WRITE ( iw, '( A, T71, I10 )' ) &
' CONTROL| Number of allocated particles ', natoms
WRITE ( iw, '( A, T71, I10 )' ) &
' CONTROL| Number of allocated particle nodes ', nnodes
WRITE ( iw, '( A, I5, T71, I10 )' ) &
' CONTROL| Number of allocated molecules on processor ', &
globenv % 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 ( 'control', 'part' )
DEALLOCATE ( struc % pnode, STAT = ios )
IF ( ios /= 0 ) CALL stop_memory ( 'control', 'pnode' )
DEALLOCATE ( struc % molecule, STAT = ios )
IF ( ios /= 0 ) CALL stop_memory ( 'control', 'molecule' )
END SUBROUTINE deallocmem
!******************************************************************************
END MODULE control_module

22
src/cp2k.F
View file

@ -27,14 +27,12 @@ PROGRAM cp2k
USE cp2k_input, ONLY : read_cp2k_section
USE environment, ONLY : initialisation, trailer
USE fist, ONLY : fist_main
USE fist_global, ONLY : set_fist_global
USE control_module, ONLY : control
USE global_types, ONLY : global_environment_type
USE library_tests,ONLY : lib_test
USE kinds, ONLY : dbl, print_kind_info
USE parallel, ONLY : start_parallel, end_parallel
USE physcon, ONLY : print_physcon
USE tbmd, ONLY : tbmd_main
USE tbmd_global, ONLY : set_tbmd_global
IMPLICIT NONE
@ -53,17 +51,19 @@ PROGRAM cp2k
IF ( globenv % ionode .AND. globenv % print_level > 4 ) &
CALL print_physcon ( globenv % scr )
IF ( globenv % program_name == 'FIST' ) THEN
CALL set_fist_global ( globenv )
CALL fist_main()
IF ( globenv % program_name == 'FIST' .OR. &
globenv % program_name == 'TBMD' ) THEN
CALL control ( globenv )
ELSE IF ( globenv % program_name == 'TBMD' ) THEN
CALL set_tbmd_global ( globenv )
CALL tbmd_main()
ELSE IF ( globenv % program_name == 'TEST' ) THEN
CALL lib_test ( globenv )
END IF
CALL trailer ( globenv )
CALL end_parallel()
CALL end_parallel ( )
END PROGRAM cp2k

106
src/dgs.F
View file

@ -6,12 +6,12 @@
MODULE dgs
USE coefficient_types, ONLY : coeff_type
USE fft_tools, ONLY : fft_radix_operations, &
FFT_RADIX_ALLOWED, FFT_RADIX_DISALLOWED, BWFFT, FWFFT, &
fft3d => fft_wrap
USE fft_tools, ONLY : fft_radix_operations, fft3d, &
FFT_RADIX_ALLOWED, FFT_RADIX_DISALLOWED, FFT_RADIX_NEXT, &
FFT_RADIX_CLOSEST, BWFFT, FWFFT
USE kinds, ONLY : dbl, sgl
USE mathconstants, ONLY : twopi, pi
USE pw_grid_types, ONLY : pw_grid_type
USE pw_grid_types, ONLY : pw_grid_type, HALFSPACE
USE pw_grids, ONLY : pw_grid_setup, pw_find_cutoff
USE pw_types, ONLY : COMPLEXDATA3D
USE pws, ONLY : fft_wrap
@ -38,46 +38,43 @@ SUBROUTINE dg_grid_setup ( box_b, npts_s, epsilon, alpha, grid_s, grid_b, &
! Arguments
INTEGER, DIMENSION ( : ), INTENT ( IN ) :: npts_s
REAL ( dbl ), INTENT ( IN ) :: epsilon, alpha
REAL ( dbl ), INTENT ( INOUT ) :: epsilon
REAL ( dbl ), INTENT ( IN ) :: alpha
TYPE ( cell_type ), INTENT ( IN ) :: box_b
CHARACTER( LEN = * ), INTENT ( IN ) :: dg_gaussian_type
TYPE ( pw_grid_type ), INTENT ( INOUT ) :: grid_s, grid_b
! Locals
INTEGER :: check
INTEGER :: nout ( 3 )
REAL ( dbl ) :: dr ( 3 ), cutoff
REAL ( dbl ) :: cell_lengths ( 3 )
TYPE ( cell_type ) :: unit_box, box_s
INTEGER :: foo ( 3 ) !*apsi
!------------------------------------------------------------------------------
CALL fft_radix_operations ( npts_s ( 1 ), check, &
operation = FFT_RADIX_ALLOWED )
IF ( check /= FFT_RADIX_ALLOWED ) THEN
CALL stop_prg ( "dg_grid_setup", "disallowed small FFT length #1" )
END IF
CALL fft_radix_operations ( npts_s ( 2 ), check, &
operation = FFT_RADIX_ALLOWED )
IF ( check /= FFT_RADIX_ALLOWED ) THEN
CALL stop_prg ( "dg_grid_setup", "disallowed small FFT length #2" )
END IF
CALL fft_radix_operations ( npts_s ( 3 ), check, &
operation = FFT_RADIX_ALLOWED )
IF ( check /= FFT_RADIX_ALLOWED ) THEN
CALL stop_prg ( "dg_grid_setup", "disallowed small FFT length #3" )
END IF
CALL fft_radix_operations ( npts_s ( 1 ), nout ( 1 ), &
operation = FFT_RADIX_NEXT )
CALL fft_radix_operations ( npts_s ( 1 ), nout ( 2 ), &
operation = FFT_RADIX_NEXT )
CALL fft_radix_operations ( npts_s ( 1 ), nout ( 3 ), &
operation = FFT_RADIX_NEXT )
CALL get_cell_param ( box_b, cell_lengths )
CALL dg_get_spacing ( npts_s, epsilon, alpha, dg_gaussian_type, dr )
CALL dg_get_spacing ( nout, epsilon, alpha, dg_gaussian_type, dr )
CALL dg_find_radix ( dr, cell_lengths, grid_b % npts )
CALL dg_get_epsilon ( nout, epsilon, alpha, dg_gaussian_type, dr )
grid_b % bounds ( 1, : ) = - grid_b % npts / 2
grid_b % bounds ( 2, : ) = + ( grid_b % npts - 1 ) / 2
grid_s % npts ( : ) = grid_s % bounds ( 2, : ) - grid_s % bounds ( 1, : ) + 1
grid_b % grid_span = HALFSPACE
grid_s % bounds ( 1, : ) = -nout ( : ) / 2
grid_s % bounds ( 2, : ) = ( +nout ( : ) - 1 ) / 2
grid_s % grid_span = HALFSPACE
grid_s % npts = nout
CALL pw_find_cutoff ( grid_b % npts, box_b, cutoff )
@ -109,8 +106,6 @@ SUBROUTINE dg_get_spacing ( npts, epsilon, alpha, dg_gaussian_type, dr )
!------------------------------------------------------------------------------
write(6,*) "dg_get_spacing, please check n/2 vs REAL(n)/2.0 !!!"
SELECT CASE ( dg_gaussian_type )
CASE ( "PME_GAUSS" ) ! use a Gaussian of form N exp(-2*alpha^2*r^2)
@ -137,18 +132,71 @@ END SUBROUTINE dg_get_spacing
!******************************************************************************
SUBROUTINE dg_get_epsilon ( npts, epsilon, alpha, dg_gaussian_type, dr )
IMPLICIT NONE
! Arguments
INTEGER, DIMENSION ( : ), INTENT ( IN ) :: npts
REAL ( dbl ), INTENT ( OUT ) :: epsilon
REAL ( dbl ), INTENT ( IN ) :: alpha
CHARACTER ( LEN = * ), INTENT ( IN ) :: dg_gaussian_type
REAL ( dbl ), DIMENSION ( : ), INTENT ( OUT ) :: dr
! Locals
REAL ( dbl ) :: alphasq, norm, rlen
!------------------------------------------------------------------------------
SELECT CASE ( dg_gaussian_type )
CASE ( "PME_GAUSS" ) ! use a Gaussian of form N exp(-2*alpha^2*r^2)
rlen = MAXVAL ( REAL ( npts, dbl ) * dr ) * 0.5_dbl
alphasq = alpha ** 2
norm = ( 2.0_dbl * alphasq / pi ) ** ( 1.5_dbl )
epsilon = norm * exp ( - 2._dbl * alphasq * rlen * rlen )
CASE ( "S" ) ! use a Gaussian of form N exp(-2*alpha^2*r^2)
rlen = MAXVAL ( REAL ( npts, dbl ) * dr ) * 0.5_dbl
alphasq = alpha ** 2
norm = ( 2.0_dbl * alphasq / pi ) ** ( 1.5_dbl )
epsilon = norm * exp ( - 2._dbl * alphasq * rlen * rlen )
CASE ( "P" )
CALL stop_prg ( "dg_get_spacing", "'P; gaussian type not defined" )
CASE DEFAULT
CALL stop_prg ( "dg_get_spacing", "no suitable gaussian type specified" )
END SELECT
END SUBROUTINE dg_get_epsilon
!******************************************************************************
SUBROUTINE dg_find_radix ( dr, cell_lengths, npts )
IMPLICIT NONE
! Arguments
REAL ( dbl ), INTENT ( IN ) :: dr ( 3 )
REAL ( dbl ), INTENT ( INOUT ) :: dr ( 3 )
REAL ( dbl ), INTENT ( IN ) :: cell_lengths ( 3 )
INTEGER, DIMENSION ( : ), INTENT ( OUT ) :: npts
! Locals
INTEGER, DIMENSION ( 3 ) :: nin
!------------------------------------------------------------------------------
npts ( : ) = NINT ( cell_lengths ( : ) / dr ( : ) )
nin ( : ) = NINT ( cell_lengths ( : ) / dr ( : ) )
CALL fft_radix_operations ( nin ( 1 ), npts ( 1 ), &
operation = FFT_RADIX_CLOSEST )
CALL fft_radix_operations ( nin ( 2 ), npts ( 2 ), &
operation = FFT_RADIX_CLOSEST )
CALL fft_radix_operations ( nin ( 3 ), npts ( 3 ), &
operation = FFT_RADIX_CLOSEST )
dr ( : ) = cell_lengths ( : ) / REAL ( npts ( : ), dbl )
END SUBROUTINE dg_find_radix

4
src/environment.F
View file

@ -15,6 +15,7 @@ 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
@ -86,6 +87,9 @@ SUBROUTINE initialisation ( globenv )
! initialize physical constants
CALL init_physcon()
! initialize FFT library
CALL init_fft ( )
END SUBROUTINE initialisation

67
src/ewalds.F
View file

@ -9,6 +9,7 @@ MODULE ewalds
USE dgs, ONLY : dg_grid_setup
USE dg_types, ONLY : dg_type
USE ewald_parameters_types, ONLY : ewald_parameters_type
USE global_types, ONLY : global_environment_type
USE kinds, ONLY : dbl
USE mathconstants, ONLY : pi, zero
USE md, ONLY : thermodynamic_type
@ -317,19 +318,18 @@ END SUBROUTINE ewald_print
!******************************************************************************
SUBROUTINE ewald_initialize ( dg, part, pnode, pnode_grp, ewald_param, box, &
thermo, iounit, ewald_grid, pme_small_grid, pme_big_grid )
thermo, ewald_grid, pme_small_grid, pme_big_grid )
IMPLICIT NONE
! Arguments
TYPE ( dg_type ), INTENT ( OUT ) :: dg
TYPE ( ewald_parameters_type ), INTENT ( IN ) :: ewald_param
TYPE ( ewald_parameters_type ), INTENT ( INOUT ) :: ewald_param
TYPE ( particle_type ), DIMENSION ( : ), INTENT ( IN ) :: part
TYPE ( particle_node_type ), DIMENSION ( : ), INTENT ( IN ) :: pnode
INTEGER, INTENT ( IN ) :: pnode_grp
TYPE ( global_environment_type ), INTENT ( IN ) :: pnode_grp
TYPE ( cell_type ), INTENT ( IN ) :: box
TYPE ( thermodynamic_type ), INTENT ( INOUT ) :: thermo
INTEGER, INTENT ( IN ) :: iounit
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
@ -342,34 +342,36 @@ SUBROUTINE ewald_initialize ( dg, part, pnode, pnode_grp, ewald_param, box, &
! parallelisation is over atoms (pnodes), so the group of processors
! has to be the same as the group for the pnodes
ewald_grp = pnode_grp
ewald_grp = pnode_grp % group
! writing output to unit iounit
iw = iounit
iw = pnode_grp % scr
natoms = SIZE ( part )
IF ( ewald_param % ewald_type /= 'NONE' ) THEN
IF ( pnode_grp % ionode ) THEN
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
WRITE ( iw, '( A,T67,A14 )' ) ' Ewald| Summation is done by:', &
ADJUSTR(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 (input) ', ewald_param % ns_max
WRITE ( iw, '( A,T71,E10.4 )' ) &
' PME| Gaussian tolerance (input) ', ewald_param % epsilon
END SELECT
ELSE
ELSE
WRITE ( iw, '( A )' ) ' No Ewald summation is performed'
WRITE ( iw, '( A, T73, A )' ) ' Ewald| ','not used'
END IF
END IF
! fire up the reciprocal space and compute self interaction and
@ -383,7 +385,8 @@ SUBROUTINE ewald_initialize ( dg, part, pnode, pnode_grp, ewald_param, box, &
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" )
IF ( pnode_grp % ionode ) &
CALL stop_prg ( "initialize_ewalds", "gmax has to be odd" )
END IF
ewald_grid % bounds ( 1, : ) = -gmax / 2
ewald_grid % bounds ( 2, : ) = +gmax / 2
@ -401,16 +404,22 @@ SUBROUTINE ewald_initialize ( dg, part, pnode, pnode_grp, ewald_param, box, &
IF ( PRESENT ( pme_small_grid ) .AND. PRESENT ( pme_big_grid ) ) THEN
npts_s ( : ) = ewald_param % ns_max
pme_small_grid % bounds ( 1, : ) = -npts_s ( : ) / 2
pme_small_grid % bounds ( 2, : ) = ( +npts_s ( : ) - 1 ) / 2
pme_small_grid % grid_span = HALFSPACE
pme_big_grid % grid_span = HALFSPACE
CALL dg_grid_setup ( box, npts_s, ewald_param % epsilon, &
ewald_param % alpha, pme_small_grid, &
pme_big_grid, ewald_param % ewald_type )
CALL pme_setup (pnode, pme_small_grid, ewald_param, dg )
IF ( pnode_grp % ionode ) THEN
WRITE ( iw, '( A,T71,E10.4 )' ) &
' PME| Gaussian tolerance (effective) ', ewald_param % epsilon
WRITE ( iw, '( A,T63,3I6 )' ) &
' PME| Small box grid ', pme_small_grid % npts
WRITE ( iw, '( A,T63,3I6 )' ) &
' PME| Full box grid ', pme_big_grid % npts
END IF
END IF
END IF
@ -420,3 +429,5 @@ END SUBROUTINE ewald_initialize
!******************************************************************************
END MODULE ewalds
!******************************************************************************

176
src/fermi.F Normal file
View file

@ -0,0 +1,176 @@
!------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!------------------------------------------------------------------------------!
!
MODULE fermi
!
!------------------------------------------------------------------------------!
USE kinds, ONLY : dbl
USE stop_program, ONLY : stop_prg
USE physcon, ONLY : boltzmann, joule
!
IMPLICIT NONE
!
PRIVATE
PUBLIC :: fermi_distribution_type, init_fermi_dist, fermi_info
!
TYPE fermi_distribution_type
REAL (dbl) :: electronic_temp
REAL (dbl) :: betael
REAL (dbl) :: mu
INTEGER :: nel, nalpha, nbeta
INTEGER :: spin_polarization
INTEGER :: multiplicity
INTEGER :: nstate, na, nb
END TYPE fermi_distribution_type
!------------------------------------------------------------------------------!
!
CONTAINS
!
!------------------------------------------------------------------------------!
SUBROUTINE init_fermi_dist(fdist,nel,etemp,lsd,mult,nstate)
IMPLICIT NONE
TYPE (fermi_distribution_type), INTENT (INOUT) :: fdist
INTEGER, OPTIONAL, INTENT (IN) :: nel
REAL (dbl), OPTIONAL, INTENT (IN) :: etemp
LOGICAL, OPTIONAL, INTENT (IN) :: lsd
INTEGER, OPTIONAL, INTENT (IN) :: mult
INTEGER, OPTIONAL, INTENT (IN) :: nstate
INTEGER :: m, isos
fdist%nel = 0
fdist%electronic_temp = 0._dbl
fdist%spin_polarization = 0
fdist%multiplicity = 1
fdist%nstate = 0
fdist%mu = 0._dbl
!..total number of electrons
IF (present(nel)) THEN
fdist%nel = nel
END IF
!..electronic temperature
IF (present(etemp)) THEN
fdist%electronic_temp = etemp
END IF
IF (fdist%electronic_temp>0._dbl) THEN
fdist%betael = joule/(boltzmann*fdist%electronic_temp)
ELSE
fdist%betael = 1.E33_dbl
END IF
!..spin polarization
IF (present(lsd)) THEN
IF (lsd) fdist%spin_polarization = 1
END IF
!..multiplicity
IF (present(mult)) THEN
IF (lsd) fdist%multiplicity = mult
END IF
!..number of electronic states
IF (present(nstate)) THEN
fdist%nstate = nstate
ELSE
fdist%nstate = 0
END IF
!..test the current setting and adjust defaults
!..number of electrons at least 1
IF (fdist%nel<1) THEN
CALL stop_prg('INIT_FERMI_DIST','Number of electrons < 1')
END IF
!..if lda only singlet states
IF (fdist%spin_polarization==0 .AND. fdist%multiplicity/=1) THEN
CALL stop_prg('INIT_FERMI_DIST','LDA: multiplicity has to be 1')
END IF
!..set number of alpha and beta electrons
IF (fdist%spin_polarization==1) THEN
m = fdist%multiplicity
IF (mod(nel+m-1,2)/=0) THEN
CALL stop_prg('INIT_FERMI_DIST','NEL incons. with multiplicity')
END IF
fdist%nalpha = (nel+m-1)/2
fdist%nbeta = nel - fdist%nalpha
ELSE
fdist%nalpha = nel
fdist%nbeta = 0
END IF
!..set number of states
IF (fdist%nstate/=0) THEN
IF (fdist%spin_polarization==1) THEN
fdist%na = fdist%nalpha
fdist%nb = fdist%nbeta
m = fdist%nstate - (fdist%na+fdist%nb)
IF (m<0) THEN
CALL stop_prg('INIT_FERMI_DIST','Not enough states')
END IF
fdist%na = fdist%na + (m+1)/2
fdist%nb = fdist%nb + m/2
ELSE
fdist%na = fdist%nstate
fdist%nb = 0
IF (fdist%nstate<(fdist%nel+1)/2) THEN
CALL stop_prg('INIT_FERMI_DIST','Not enough states')
END IF
END IF
ELSE
IF (fdist%electronic_temp==0._dbl) THEN
IF (fdist%spin_polarization==1) THEN
fdist%na = fdist%nalpha
fdist%nb = fdist%nbeta
fdist%nstate = fdist%na + fdist%nb
ELSE
fdist%nstate = (fdist%nel+1)/2
fdist%na = fdist%nstate
fdist%nb = 0
END IF
ELSE
CALL stop_prg('INIT_FERMI_DIST','If etemp /= 0', &
' number of states has to be specified')
END IF
END IF
END SUBROUTINE init_fermi_dist
!------------------------------------------------------------------------------!
SUBROUTINE fermi_info(fdist,plevel,punit)
IMPLICIT NONE
TYPE (fermi_distribution_type), INTENT (IN) :: fdist
INTEGER :: plevel, punit
CHARACTER (len=10) :: mul(1:8)
mul = ' '
mul(1) = 'singlet'
mul(2) = 'doublet'
mul(3) = 'triplet'
mul(4) = 'quartet'
mul(5) = 'quintet'
mul(6) = 'sextet'
mul(7) = 'septet'
mul(8) = 'octet'
IF (plevel>0) THEN
WRITE (punit,*)
WRITE (punit,'(A,T71,I10)') ' FERMI| Total number of electrons ', &
fdist%nel
IF (fdist%spin_polarization==1) THEN
WRITE (punit,'(A,T71,I10)') ' FERMI| Number of alpha electrons ', &
fdist%nalpha
WRITE (punit,'(A,T71,I10)') ' FERMI| Number of beta electrons ', &
fdist%nbeta
END IF
WRITE (punit,'(A,T71,A)') ' FERMI| Multiplicity', &
adjustr(mul(fdist%multiplicity))
WRITE (punit,'(A,T71,F10.2)') ' FERMI| Electronic temperature [K]', &
fdist%electronic_temp
WRITE (punit,'(A,T71,I10)') &
' FERMI| Total number of electronic states ', fdist%nstate
IF (fdist%spin_polarization==1) THEN
WRITE (punit,'(A,T71,I10)') &
' FERMI| Number of alpha electronic states ', fdist%na
WRITE (punit,'(A,T71,I10)') &
' FERMI| Number of beta electronic states ', fdist%nb
END IF
END IF
END SUBROUTINE fermi_info
!------------------------------------------------------------------------------!
END MODULE fermi
!------------------------------------------------------------------------------!

992
src/fft_tools.F
View file

File diff suppressed because it is too large Load diff

155
src/fftsg_lib.F Normal file
View file

@ -0,0 +1,155 @@
!-----------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
MODULE fftsg_lib
USE kinds, ONLY: dbl
USE stop_program, ONLY : stop_memory, stop_prg
IMPLICIT NONE
PRIVATE
PUBLIC :: fft3d
PUBLIC :: fft_get_lengths
!******************************************************************************
CONTAINS
!******************************************************************************
!! Give the allowed lengths of FFT's '''
SUBROUTINE fft_get_lengths ( data, max_length )
IMPLICIT NONE
! Arguments
INTEGER, INTENT ( IN ) :: max_length
INTEGER, DIMENSION ( : ), POINTER :: data
! Locals
INTEGER, PARAMETER :: rlen = 82
INTEGER, DIMENSION ( rlen ), PARAMETER :: radix = &
(/ 2, 4, 5, 6, 8, 9, 12, 15, 16, 18, 20, 24, 25, 27, 30, 32, 36, 40, &
45, 48, 54, 60, 64, 72, 75, 80, 81, 90, 96, 100, 108, 120, 125, 128, &
135, 144, 150, 160, 162, 180, 192, 200, 216, 225, 240, 243, 256, 270, &
288, 300, 320, 324, 360, 375, 384, 400, 405, 432, 450, 480, 486, 500, &
512, 540, 576, 600, 625, 640, 648, 675, 720, 729, 750, 768, 800, 810, &
864, 900, 960, 972, 1000, 1024 /)
INTEGER :: i, allocstat, ndata
!------------------------------------------------------------------------------
ndata = 0
DO i = 1, rlen
IF ( radix ( i ) > max_length ) EXIT
ndata = ndata + 1
END DO
ALLOCATE ( data ( ndata ), STAT = allocstat )
IF ( allocstat /= 0 ) THEN
CALL stop_memory ( "fft_get_lengths", "data", ndata )
END IF
data ( 1:ndata ) = radix ( 1:ndata )
END SUBROUTINE fft_get_lengths
!******************************************************************************
! routine with wrapper for all fft calls:
! Does transform with exp(+ig.r*sign):
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
COMPLEX ( dbl ), DIMENSION(:), ALLOCATABLE :: xf, yf
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 defined ( __FFTSG )
IF ( fft_in_place ) THEN
ALLOCATE ( xf ( ldx*ldy*ldz ), STAT = ierr )
IF ( ierr /= 0 ) call stop_memory ( "fft3d", "xf", ldx*ldy*ldz )
ALLOCATE ( yf ( ldx*ldy*ldz ), STAT = ierr )
IF ( ierr /= 0 ) call stop_memory ( "fft3d", "yf", ldx*ldy*ldz )
CALL mltfftsg ( 'N', 'T', zg, ldx, ldy*ldz, xf, ldy*ldz, ldx, nx, &
ldy*ldz, sign_fft, 1._dbl )
CALL mltfftsg ( 'N', 'T', xf, ldy, ldx*ldz, yf, ldx*ldz, ldy, ny, &
ldx*ldz, sign_fft, 1._dbl )
CALL mltfftsg ( 'N', 'T', yf, ldz, ldy*ldx, zg, ldy*ldx, ldz, nz, &
ldy*ldx, sign_fft, scale)
DEALLOCATE ( xf, STAT = ierr )
IF ( ierr /= 0 ) call stop_memory ( "fft3d", "xf" )
DEALLOCATE ( yf, STAT = ierr )
IF ( ierr /= 0 ) call stop_memory ( "fft3d", "yf" )
ELSE
ldox = SIZE ( zg_out (:,1,1) )
ldoy = SIZE ( zg_out (1,:,1) )
ldoz = SIZE ( zg_out (1,1,:) )
ALLOCATE ( xf ( ldx*ldy*ldz ), STAT = ierr )
IF ( ierr /= 0 ) call stop_memory ( "fft3d", "xf", ldx*ldy*ldz )
CALL mltfftsg ( 'N', 'T', zg, ldx, ldy*ldz, zg_out, ldy*ldz, ldx, nx, &
ldy*ldz, sign_fft, 1._dbl )
CALL mltfftsg ( 'N', 'T', zg_out, ldy, ldx*ldz, xf, ldx*ldz, ldy, ny, &
ldx*ldz, sign_fft, 1._dbl )
CALL mltfftsg ( 'N', 'T', xf, ldz, ldy*ldx, zg_out, ldy*ldx, ldz, nz, &
ldy*ldx, sign_fft, scale)
DEALLOCATE ( xf, STAT = ierr )
IF ( ierr /= 0 ) call stop_memory ( "fft3d", "xf" )
END IF
#else
fsign = 0
#endif
END SUBROUTINE fft3d
!******************************************************************************
END MODULE fftsg_lib
!******************************************************************************

313
src/fftw_lib.F Normal file
View file

@ -0,0 +1,313 @@
!-----------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
MODULE fftw_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
!*apsi 220500 This is actually the file fftw_f77.i ...
! This file contains PARAMETER statements for various constants
! that can be passed to FFTW routines. You should include
! this file in any FORTRAN program that calls the fftw_f77
! routines (either directly or with an #include statement
! if you use the C preprocessor).
integer FFTW_FORWARD,FFTW_BACKWARD
parameter (FFTW_FORWARD=-1,FFTW_BACKWARD=1)
integer FFTW_REAL_TO_COMPLEX,FFTW_COMPLEX_TO_REAL
parameter (FFTW_REAL_TO_COMPLEX=-1,FFTW_COMPLEX_TO_REAL=1)
integer FFTW_ESTIMATE,FFTW_MEASURE
parameter (FFTW_ESTIMATE=0,FFTW_MEASURE=1)
integer FFTW_OUT_OF_PLACE,FFTW_IN_PLACE,FFTW_USE_WISDOM
parameter (FFTW_OUT_OF_PLACE=0)
parameter (FFTW_IN_PLACE=8,FFTW_USE_WISDOM=16)
integer FFTW_THREADSAFE
parameter (FFTW_THREADSAFE=128)
! Constants for the MPI wrappers:
integer FFTW_TRANSPOSED_ORDER, FFTW_NORMAL_ORDER
integer FFTW_SCRAMBLED_INPUT, FFTW_SCRAMBLED_OUTPUT
parameter(FFTW_TRANSPOSED_ORDER=1, FFTW_NORMAL_ORDER=0)
parameter(FFTW_SCRAMBLED_INPUT=8192)
parameter(FFTW_SCRAMBLED_OUTPUT=16384)
!******************************************************************************
CONTAINS
!******************************************************************************
!! Give the allowed lengths of FFT's '''
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 = 3
maxn_fives = 2
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
!******************************************************************************
! routine with wrapper for all fft calls:
! Does transform with exp(+ig.r*sign):
SUBROUTINE fft3d ( fsign, scale, n, zg, zg_out )
IMPLICIT NONE
! Arguments
INTEGER, INTENT ( INOUT ) :: fsign
REAL ( dbl ), INTENT ( IN ) :: scale
INTEGER, DIMENSION ( : ), INTENT ( IN ) :: n
COMPLEX ( dbl ), DIMENSION(:,:,:), INTENT ( INOUT ) :: zg
COMPLEX ( dbl ), DIMENSION(:,:,:), INTENT ( INOUT ), OPTIONAL, TARGET :: zg_out
! Locals
INTEGER, SAVE :: n1a_save = -1, n2a_save = -1, n3a_save = -1
INTEGER, SAVE :: n1b_save = -1, n2b_save = -1, n3b_save = -1
INTEGER, SAVE :: n1c_save = -1, n2c_save = -1, n3c_save = -1
LOGICAL, SAVE :: ffta_in_place = .TRUE., fftb_in_place = .TRUE.
LOGICAL, SAVE :: fftc_in_place = .TRUE.
LOGICAL :: fft_in_place
INTEGER :: sign_fft, n1, n2, n3
INTEGER ( KIND = 8 ), SAVE :: plan_a_fw, plan_a_bw
INTEGER ( KIND = 8 ), SAVE :: plan_b_fw, plan_b_bw
INTEGER ( KIND = 8 ), SAVE :: plan_c_fw, plan_c_bw
REAL ( dbl ) :: norm
! Just due to DEC... apsi
COMPLEX ( dbl ), DIMENSION(:,:,:), POINTER :: zgout
COMPLEX ( dbl ), DIMENSION(1,1,1), TARGET :: zgout_tmp
!------------------------------------------------------------------------------
norm = scale
n1 = n(1)
n2 = n(2)
n3 = n(3)
IF ( PRESENT ( zg_out ) ) THEN
fft_in_place = .false.
zgout => zg_out
ELSE
fft_in_place = .true.
zgout => zgout_tmp
END IF
sign_fft = fsign
#if defined ( __FFTW )
IF ( n1a_save == n1 .AND. n2a_save == n2 .AND. n3a_save == n3 &
.AND. ( fft_in_place .eqv. ffta_in_place ) ) THEN
IF ( sign_fft == +1 ) THEN
CALL fftwnd_f77_one ( plan_a_fw, zg, zgout )
ELSE
CALL fftwnd_f77_one ( plan_a_bw, zg, zgout )
END IF
ELSE IF ( n1b_save == n1 .AND. n2b_save == n2 .AND. n3b_save == n3 &
.AND. ( fft_in_place .eqv. fftb_in_place ) ) THEN
IF ( sign_fft == +1 ) THEN
CALL fftwnd_f77_one ( plan_b_fw, zg, zgout )
ELSE
CALL fftwnd_f77_one ( plan_b_bw, zg, zgout )
END IF
ELSE IF ( n1c_save == n1 .AND. n2c_save == n2 .AND. n3c_save == n3 &
.AND. ( fft_in_place .eqv. fftc_in_place ) ) THEN
IF ( sign_fft == +1 ) THEN
CALL fftwnd_f77_one ( plan_c_fw, zg, zgout )
ELSE
CALL fftwnd_f77_one ( plan_c_bw, zg, zgout )
END IF
ELSE IF ( n1a_save == -1 .OR. &
( n1a_save == n1 .AND. n2a_save == n2 .AND. n3a_save == n3 &
.AND. ( fft_in_place .neqv. ffta_in_place ) .AND. &
( n1b_save /= -1 .AND. n1c_save /= -1 ) ) ) THEN ! Initialise 'a'
IF ( fft_in_place ) THEN
CALL fftw3d_f77_create_plan ( plan_a_fw, n1, n2, n3, FFTW_FORWARD, &
FFTW_ESTIMATE + FFTW_IN_PLACE )
CALL fftw3d_f77_create_plan ( plan_a_bw, n1, n2, n3, FFTW_BACKWARD, &
FFTW_ESTIMATE + FFTW_IN_PLACE )
ELSE
CALL fftw3d_f77_create_plan ( plan_a_fw, n1, n2, n3, FFTW_FORWARD, &
FFTW_ESTIMATE + FFTW_OUT_OF_PLACE )
CALL fftw3d_f77_create_plan ( plan_a_bw, n1, n2, n3, FFTW_BACKWARD, &
FFTW_ESTIMATE + FFTW_OUT_OF_PLACE )
END IF
n1a_save = n1
n2a_save = n2
n3a_save = n3
ffta_in_place = fft_in_place
IF ( sign_fft == +1 ) THEN
CALL fftwnd_f77_one ( plan_a_fw, zg, zgout )
ELSE
CALL fftwnd_f77_one ( plan_a_bw, zg, zgout )
END IF
ELSE IF ( n1b_save == -1 .OR. &
( n1b_save == n1 .AND. n2b_save == n2 .AND. n3b_save == n3 &
.AND. ( fft_in_place .neqv. fftb_in_place ) .AND. &
n1c_save /= -1 ) ) THEN ! Initialise 'a'
IF ( fft_in_place ) THEN
CALL fftw3d_f77_create_plan ( plan_b_fw, n1, n2, n3, FFTW_FORWARD, &
FFTW_ESTIMATE + FFTW_IN_PLACE )
CALL fftw3d_f77_create_plan ( plan_b_bw, n1, n2, n3, FFTW_BACKWARD, &
FFTW_ESTIMATE + FFTW_IN_PLACE )
ELSE
CALL fftw3d_f77_create_plan ( plan_b_fw, n1, n2, n3, FFTW_FORWARD, &
FFTW_ESTIMATE + FFTW_OUT_OF_PLACE )
CALL fftw3d_f77_create_plan ( plan_b_bw, n1, n2, n3, FFTW_BACKWARD, &
FFTW_ESTIMATE + FFTW_OUT_OF_PLACE )
END IF
n1b_save = n1
n2b_save = n2
n3b_save = n3
fftb_in_place = fft_in_place
IF ( sign_fft == +1 ) THEN
CALL fftwnd_f77_one ( plan_b_fw, zg, zgout )
ELSE
CALL fftwnd_f77_one ( plan_b_bw, zg, zgout )
END IF
ELSE ! Initialise 'c'
IF ( fft_in_place ) THEN
CALL fftw3d_f77_create_plan ( plan_c_fw, n1, n2, n3, FFTW_FORWARD, &
FFTW_ESTIMATE + FFTW_IN_PLACE )
CALL fftw3d_f77_create_plan ( plan_c_bw, n1, n2, n3, FFTW_BACKWARD, &
FFTW_ESTIMATE + FFTW_IN_PLACE )
ELSE
CALL fftw3d_f77_create_plan ( plan_c_fw, n1, n2, n3, FFTW_FORWARD, &
FFTW_ESTIMATE + FFTW_OUT_OF_PLACE )
CALL fftw3d_f77_create_plan ( plan_c_bw, n1, n2, n3, FFTW_BACKWARD, &
FFTW_ESTIMATE + FFTW_OUT_OF_PLACE )
END IF
n1c_save = n1
n2c_save = n2
n3c_save = n3
fftc_in_place = fft_in_place
IF ( sign_fft == +1 ) THEN
CALL fftwnd_f77_one ( plan_c_fw, zg, zgout )
ELSE
CALL fftwnd_f77_one ( plan_c_bw, zg, zgout )
END IF
END IF
#else
fsign = 0
#endif
IF ( norm /= 1._dbl ) THEN
IF ( fft_in_place ) THEN
zg = zg * norm
ELSE
zgout = zgout * norm
END IF
END IF
END SUBROUTINE fft3d
!******************************************************************************
END MODULE fftw_lib
!******************************************************************************

225
src/fist_debug.F
View file

@ -6,23 +6,23 @@
MODULE fist_debug
USE ewald_parameters_types, ONLY : ewald_parameters_type
USE global_types, ONLY : global_environment_type
USE kinds, ONLY : dbl
USE md, ONLY : thermodynamic_type
USE molecule_types, ONLY : molecule_structure_type, particle_node_type
USE particle_types, ONLY : particle_type
USE pair_potential, ONLY : potentialparm_type, ener_coul
USE fist_force, ONLY : force_control, debug_variables_type
USE fist_force_numer, ONLY : force_bond_numer, force_bend_numer, &
force_nonbond_numer, force_recip_numer, pvbond_numer, pvbend_numer, &
ptens_numer, pvg_numer, potential_g_numer, de_g_numer, &
energy_recip_numer
USE fist_force, ONLY : force_control, debug_variables_type
USE global_types, ONLY : global_environment_type
USE kinds, ONLY : dbl
USE linklists, ONLY : bonds, bends, torsions, dist_constraints, &
g3x3_constraints
USE pw_grid_types, ONLY : pw_grid_type
USE pw_grids, ONLY : pw_grid_setup
USE md, ONLY : thermodynamic_type
USE molecule_types, ONLY : molecule_structure_type, particle_node_type
USE particle_types, ONLY : particle_type
USE pair_potential, ONLY : potentialparm_type, ener_coul
USE pw_grid_types, ONLY : pw_grid_type, HALFSPACE
USE pw_grids, ONLY : pw_find_cutoff, pw_grid_setup
USE simulation_cell, ONLY : cell_type
USE stop_program, ONLY : stop_memory
USE stop_program, ONLY : stop_memory, stop_prg
IMPLICIT NONE
@ -53,16 +53,16 @@ SUBROUTINE debug_control ( globenv, ewald_param, part, pnode, molecule, box, &
! Locals
TYPE ( debug_variables_type ) :: dbg
TYPE ( pw_grid_type ) :: pw_grid
REAL ( dbl ), DIMENSION ( :, : ), ALLOCATABLE :: f_numer
REAL ( dbl ), DIMENSION ( :, : ), ALLOCATABLE :: rel, diff
INTEGER :: iflag, i, natoms, isos, iatom, iw, ir
TYPE ( pw_grid_type ) :: ewald_grid
INTEGER :: iflag, i, natoms, isos, iatom, iw, ir, npts_s(3), gmax
INTEGER, DIMENSION ( 2 ) :: dum
REAL ( dbl ) :: delta, energy_numer
REAL ( dbl ) :: delta, energy_numer, cutoff
REAL ( dbl ) :: e_numer, pv_test ( 3, 3 )
REAL ( dbl ) :: err1, numer, denom1, vec ( 3 ), e_bc, e_real, energy_tot
REAL ( dbl ) :: denom2, err2
REAL ( dbl ), DIMENSION ( :, : ), ALLOCATABLE :: f_bc, f_real
REAL ( dbl ), DIMENSION ( :, : ), ALLOCATABLE :: f_numer
REAL ( dbl ), DIMENSION ( :, : ), ALLOCATABLE :: rel, diff
!------------------------------------------------------------------------------
@ -113,8 +113,8 @@ SUBROUTINE debug_control ( globenv, ewald_param, part, pnode, molecule, box, &
ELSE
WRITE ( iw, '( A, T71, F10.6 )' ) ' DELTA = ', delta
END IF
CALL force_nonbond_numer ( pnode, box, potparm, delta, f_numer, &
energy_numer )
CALL force_nonbond_numer ( ewald_param, pnode, box, potparm, &
delta, f_numer, energy_numer )
WRITE ( iw, '( A, T61, E20.14 )' ) ' NON BOND NUMER ENERGY = ', &
energy_numer
WRITE ( iw, '( A, T61, E20.14 )' ) ' NON BOND ANAL ENERGY = ', &
@ -314,101 +314,112 @@ SUBROUTINE debug_control ( globenv, ewald_param, part, pnode, molecule, box, &
END IF
! Debug g-space
! IF ( ewald_param % ewald_type /= 'NONE' ) THEN
! WRITE ( iw, '( A )' ) ' DO YOU WANT TO DEBUG YOUR G-SPACE (1=yes)?'
! READ ( ir, * ) iflag
! IF ( iflag == 1 ) THEN
! CALL pw_grid_setup ( box, pw_grid, ewald_param % gmax )
! WRITE ( iw, '( A )' ) &
! ' 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= ', &
! pw_grid % ngpts_cut
! WRITE ( iw, '( A, T71, F10.4 )' ) 'ALPHA= ', &
! ewald_param % alpha
! CALL energy_recip_numer ( pnode, box, potparm, pw_grid, &
! energy_numer, pw_grid % ngpts_cut )
! WRITE ( iw, '( A, T61, G20.14 )' ) 'G-SPACE ANAL ENERGY = ', &
! dbg % pot_g
! WRITE ( iw, '( A, T61, G20.14 )' ) &
! 'G-SPACE NUMERICAL ENERGY = ', energy_numer
! END IF
! WRITE ( iw, '( A )' ) &
! ' DO YOU WANT TO DEBUG YOUR G-SPACE FORCES (1=yes)?'
! READ ( ir, * ) iflag
! IF ( iflag == 1 ) THEN
! 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
! WRITE ( iw, '( A, T71, F10.6 )' ) &
! ' DELTA (changed to default) = ', delta
! ELSE
! WRITE ( iw, '( A, T71, F10.6 )' ) ' DELTA = ', delta
! END IF
!
! CALL force_recip_numer ( pnode, box, potparm, gvec, delta, f_numer )
!
IF ( ewald_param % ewald_type /= 'NONE' ) THEN
WRITE ( iw, '( A )' ) ' DO YOU WANT TO DEBUG YOUR G-SPACE (1=yes)?'
READ ( ir, * ) iflag
IF ( iflag == 1 ) THEN
gmax = ewald_param % gmax
ewald_grid % bounds ( 1, : ) = -gmax / 2
ewald_grid % bounds ( 2, : ) = +gmax / 2
npts_s = (/ gmax, gmax, gmax /)
ewald_grid % grid_span = HALFSPACE
CALL pw_find_cutoff ( npts_s, box, cutoff )
CALL pw_grid_setup( box, ewald_grid, cutoff)
WRITE ( iw, '( A )' ) &
' 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= ', &
ewald_grid % ngpts_cut
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 = ', &
dbg % pot_g
WRITE ( iw, '( A, T61, G20.14 )' ) &
'G-SPACE NUMERICAL ENERGY = ', energy_numer
END IF
WRITE ( iw, '( A )' ) &
' DO YOU WANT TO DEBUG YOUR G-SPACE FORCES (1=yes)?'
READ ( ir, * ) iflag
IF ( iflag == 1 ) THEN
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
WRITE ( iw, '( A, T71, F10.6 )' ) &
' DELTA (changed to default) = ', delta
ELSE
WRITE ( iw, '( A, T71, F10.6 )' ) ' DELTA = ', delta
END IF
CALL force_recip_numer ( ewald_param, pnode, box, ewald_grid, &
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
!
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 )' ) &
! 'MAXIMUM ABSOLUTE ERROR = ', maxval(diff)
! WRITE ( iw, '( A,T61,G20.14 )' ) &
! 'MAXIMUM RELATIVE ERROR = ', maxval(rel)
WRITE ( iw, '( A,T61,G20.14 )' ) &
'MAXIMUM ABSOLUTE ERROR = ', maxval(diff)
WRITE ( iw, '( A,T61,G20.14 )' ) &
'MAXIMUM RELATIVE ERROR = ', maxval(rel)
!
! write out the particle number and forces of
! the max absolute and relative error
! 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 =', &
! dbg % f_g(:,dum(2))
! WRITE ( iw, '( A,T21,3G20.14 )' ) ' 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 =', &
! dbg % f_g(:,dum(2))
! WRITE ( iw, '( A,T21,3G20.14 )' ) ' F_NUMR G-SPACE =', &
! f_numer(:,dum(2))
! END IF
!
! WRITE ( iw, '( A )' ) &
! ' DO YOU WANT TO DEBUG YOUR G-SPACE VIRIAL (1=yes)?'
! READ ( ir, * ) iflag
! IF ( iflag == 1 ) THEN
! write out the particle number and forces of the max absolute
! and relative error
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 =', &
dbg % f_g(:,dum(2))
WRITE ( iw, '( A,T21,3G20.14 )' ) ' 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 =', &
dbg % f_g(:,dum(2))
WRITE ( iw, '( A,T21,3G20.14 )' ) ' F_NUMR G-SPACE =', &
f_numer(:,dum(2))
END IF
WRITE ( iw, '( A )' ) &
' DO YOU WANT TO DEBUG YOUR G-SPACE VIRIAL (1=yes)?'
READ ( ir, * ) iflag
IF ( iflag == 1 ) THEN
!
! get numerical virial
! 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
! WRITE ( iw, '( A,T71,F10.6 )' ) &
! ' DELTA (changed to default) = ', delta
! ELSE
! WRITE ( iw, '( A,T71,F10.6 )' ) ' DELTA = ', delta
! END IF
! CALL pvg_numer ( pnode, box, potparm, gvec, pv_test, delta )
!
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
WRITE ( iw, '( A,T71,F10.6 )' ) &
' DELTA (changed to default) = ', delta
ELSE
WRITE ( iw, '( A,T71,F10.6 )' ) ' DELTA = ', delta
END IF
CALL pvg_numer ( ewald_param, pnode, box, &
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,:)
! END DO
! WRITE ( iw, '( A )' ) ' PV G-SPACE'
! DO i = 1, 3
! WRITE ( iw, '( T21,3G20.14 )' ) dbg % pv_g(i,:)
! END DO
! END IF
! END IF
! END IF
!
WRITE ( iw, '( A )' ) ' PV G-SPACE NUMERICAL'
DO i = 1, 3
WRITE ( iw, '( T21,3G20.14 )' ) pv_test(i,:)
END DO
WRITE ( iw, '( A )' ) ' PV G-SPACE'
DO i = 1, 3
WRITE ( iw, '( T21,3G20.14 )' ) dbg % pv_g(i,:)
END DO
END IF
END IF
END IF
! Debug real-space virial
WRITE ( iw, '( A )' ) &
' DO YOU WANT TO DEBUG YOUR NON BOND (REAL SPACE) VIRIAL (1=yes)?'
@ -425,7 +436,7 @@ SUBROUTINE debug_control ( globenv, ewald_param, part, pnode, molecule, box, &
ELSE
WRITE ( iw, '( A, T71, F10.6 )' ) ' DELTA = ', delta
END IF
CALL ptens_numer ( pnode, box, potparm, pv_test, delta )
CALL ptens_numer ( ewald_param, pnode, box, potparm, pv_test, delta )
!
! writing out the virials
WRITE ( iw, '( A )' ) ' PV NUMERICAL'

52
src/fist_force.F
View file

@ -61,7 +61,7 @@ SUBROUTINE force_control ( molecule, pnode, part, box, thermo, &
TYPE ( debug_variables_type ), INTENT ( OUT ), OPTIONAL :: debug
! Locals
INTEGER :: id, i, natoms, nnodes, handle, isos
INTEGER :: id, i, ii, natoms, nnodes, handle, isos
REAL ( dbl ) :: pot_nonbond, pot_bond, pot_bend, vg_coulomb
REAL ( dbl ), DIMENSION ( :,: ), ALLOCATABLE, SAVE :: f_nonbond
REAL ( dbl ), DIMENSION ( 3,3 ) :: pv_nonbond, pv_bond, pv_bend
@ -80,21 +80,27 @@ SUBROUTINE force_control ( molecule, pnode, part, box, thermo, &
natoms = SIZE ( part )
isos = 0
first_time = .NOT. ALLOCATED ( f_nonbond )
IF ( .NOT. ALLOCATED ( f_nonbond ) ) &
ALLOCATE ( f_nonbond ( 3,natoms ), STAT = isos )
IF ( isos /= 0 ) &
IF ( .NOT. ALLOCATED ( f_nonbond ) ) THEN
ALLOCATE ( f_nonbond ( 3,natoms ), STAT = isos )
IF ( isos /= 0 ) &
CALL stop_memory ( 'force_control', 'f_nonbond', 3 * natoms )
ELSE IF ( SIZE ( f_nonbond ( 1, : ) ) < natoms ) THEN
DEALLOCATE ( f_nonbond, STAT = isos )
IF ( isos /= 0 ) CALL stop_memory ( 'force_control', 'f_nonbond' )
ALLOCATE ( f_nonbond ( 3,natoms ), STAT = isos )
IF ( isos /= 0 ) &
CALL stop_memory ( 'force_control', 'f_nonbond', 3 * natoms )
END IF
! initialize ewalds
IF ( first_time ) THEN
SELECT CASE ( ewald_param % ewald_type )
CASE ( 'EWALD_GAUSS' )
CALL ewald_initialize ( dg, part, pnode, fc_global % group, &
ewald_param, box, thermo, fc_global % scr, &
ewald_grid = grid_ewald )
CALL ewald_initialize ( dg, part, pnode, fc_global, &
ewald_param, box, thermo, ewald_grid = grid_ewald )
CASE ( 'PME_GAUSS' )
CALL ewald_initialize ( dg, part, pnode, fc_global % group, &
ewald_param, box, thermo, fc_global % scr, &
CALL ewald_initialize ( dg, part, pnode, fc_global, &
ewald_param, box, thermo, &
pme_small_grid = grid_s, pme_big_grid = grid_b )
END SELECT
END IF
@ -200,19 +206,19 @@ SUBROUTINE force_control ( molecule, pnode, part, box, thermo, &
! add up all the forces
! nonbonded forces might be claculated for atoms not on this node
! ewald forces are strictly local -> sum only over pnode
! We first sum the forces in f_nonbond, this allows for a more efficient
! global sum in the parallel code and in the end copy them back to part
DO i = 1, natoms
part ( i ) % f ( 1 ) = part ( i ) % f ( 1 ) + f_nonbond ( 1, i )
part ( i ) % f ( 2 ) = part ( i ) % f ( 2 ) + f_nonbond ( 2, i )
part ( i ) % f ( 3 ) = part ( i ) % f ( 3 ) + f_nonbond ( 3, i )
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 )
END DO
IF ( ewald_param % ewald_type /= 'NONE' ) THEN
DO i = 1, nnodes
pnode ( i ) % p % f ( 1 ) = pnode ( i ) % p % f ( 1 ) &
+ fg_coulomb ( 1, i )
pnode ( i ) % p % f ( 2 ) = pnode ( i ) % p % f ( 2 ) &
+ fg_coulomb ( 2, i )
pnode ( i ) % p % f ( 3 ) = pnode ( i ) % p % f ( 3 ) &
+ fg_coulomb ( 3, i )
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 )
END DO
END IF
@ -265,10 +271,14 @@ SUBROUTINE force_control ( molecule, pnode, part, box, thermo, &
IF ( isos /= 0 ) CALL stop_memory ( 'force_control', 'fg_coulomb' )
#if defined ( __parallel )
DO i = 1, natoms
CALL mp_sum ( part ( i ) % f, fc_global % group )
END DO
CALL mp_sum ( f_nonbond, fc_global % group )
#endif
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 )

293
src/fist_force_numer.F
View file

@ -14,6 +14,7 @@ MODULE fist_force_numer
USE mol_force, ONLY : force_bonds, force_bends
USE pair_potential, ONLY : potential_f, potentialparm_type
USE pw_grid_types, ONLY : pw_grid_type
USE pw_grids, ONLY : pw_find_cutoff, pw_grid_setup
USE simulation_cell, ONLY : cell_type, pbc, get_hinv, get_cell_param
USE stop_program, ONLY : stop_memory
@ -197,8 +198,8 @@ END SUBROUTINE force_bend_numer
!******************************************************************************
SUBROUTINE force_nonbond_numer ( pnode, box, potparm, numerical_shift, &
f_numer, energy_numer )
SUBROUTINE force_nonbond_numer ( ewald_param, pnode, box, potparm, &
numerical_shift, f_numer, energy_numer )
!
! Calculates the force and the potential of the minimum image, and
@ -207,6 +208,7 @@ SUBROUTINE force_nonbond_numer ( pnode, box, potparm, numerical_shift, &
IMPLICIT NONE
! Arguments
TYPE ( ewald_parameters_type ), INTENT ( IN ) :: ewald_param
TYPE ( particle_node_type ), DIMENSION ( : ), INTENT ( IN ) :: pnode
TYPE ( cell_type ), INTENT ( INOUT ) :: box
TYPE ( potentialparm_type ), DIMENSION ( :, : ), INTENT ( IN ) :: potparm
@ -298,17 +300,29 @@ SUBROUTINE force_nonbond_numer ( pnode, box, potparm, numerical_shift, &
delta(id) = numerical_shift
rij_minus_delta_sq = dot_product(rij-delta,rij-delta)
rij_plus_delta_sq = dot_product(rij+delta,rij+delta)
CALL potential_f ( rij_minus_delta_sq, potparm, qi, qj, &
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, &
CALL potential_f(rij_plus_delta_sq,potparm,qi,qj, &
iatomtype, jatomtype, energy_plus )
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
f_numer(id,i) = f_numer(id,i) + (energy_plus-energy_minus)
f_numer(id,j) = f_numer(id,j) - (energy_plus-energy_minus)
delta = 0.0_dbl
END DO DIM_LOOP
CALL potential_f ( rijsq, potparm, qi, qj, &
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
energy_numer = energy_numer + energy
END IF NOTMATCH
END IF
@ -339,10 +353,17 @@ SUBROUTINE force_nonbond_numer ( pnode, box, potparm, numerical_shift, &
delta(id) = numerical_shift
rij_minus_delta_sq = dot_product(rij-delta,rij-delta)
rij_plus_delta_sq = dot_product(rij+delta,rij+delta)
CALL potential_f(rij_minus_delta_sq,potparm,qi,qj, &
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, &
CALL potential_f(rij_plus_delta_sq,potparm,qi,qj, &
iatomtype, jatomtype, energy_plus )
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
f_numer(id,i) = f_numer(id,i) &
+ (energy_plus-energy_minus)
f_numer(id,j) = f_numer(id,j) &
@ -350,8 +371,13 @@ SUBROUTINE force_nonbond_numer ( pnode, box, potparm, numerical_shift, &
delta = 0.0_dbl
END DO DIM_LOOP2
CALL potential_f ( rijsq, potparm, qi, qj, &
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
energy_numer = energy_numer + energy
END IF
END DO
@ -377,10 +403,17 @@ SUBROUTINE force_nonbond_numer ( pnode, box, potparm, numerical_shift, &
delta(id) = numerical_shift
rij_minus_delta_sq = dot_product(rij-delta,rij-delta)
rij_plus_delta_sq = dot_product(rij+delta,rij+delta)
CALL potential_f(rij_minus_delta_sq,potparm,qi,qi, &
iatomtype,iatomtype,energy_minus)
CALL potential_f(rij_plus_delta_sq,potparm,qi,qi, &
iatomtype,iatomtype,energy_plus)
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
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
f_numer(id,i) = f_numer(id,i) &
+ (energy_plus-energy_minus)
f_numer(id,i) = f_numer(id,i) &
@ -388,8 +421,13 @@ SUBROUTINE force_nonbond_numer ( pnode, box, potparm, numerical_shift, &
delta = 0.0_dbl
END DO DIM_LOOP3
CALL potential_f ( rijsq, potparm, qi, qi, &
iatomtype,iatomtype, energy )
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
energy_numer = energy_numer + energy
END IF
END DO
@ -411,16 +449,15 @@ END SUBROUTINE force_nonbond_numer
!******************************************************************************
SUBROUTINE force_recip_numer ( ewald_param, pnode, box, potparm, &
pw_grid, numerical_shift, f_numer )
SUBROUTINE force_recip_numer ( ewald_param, pnode, box, pw_grid, &
numerical_shift, f_numer )
IMPLICIT NONE
! Arguments
TYPE ( cell_type ), INTENT ( IN ) :: box
TYPE ( ewald_parameters_type ), INTENT ( IN ) :: ewald_param
TYPE ( particle_node_type ), DIMENSION ( : ), INTENT ( IN ) :: pnode
TYPE ( cell_type ), INTENT ( IN ) :: box
TYPE ( potentialparm_type ), DIMENSION ( :, : ), INTENT ( IN ) :: potparm
TYPE ( pw_grid_type ), INTENT ( IN ) :: pw_grid
REAL ( dbl ), INTENT ( IN ) :: numerical_shift
REAL ( dbl ), DIMENSION ( :, : ), INTENT ( INOUT ) :: f_numer
@ -428,48 +465,30 @@ SUBROUTINE force_recip_numer ( ewald_param, pnode, box, potparm, &
! Locals
COMPLEX ( dbl ), ALLOCATABLE, DIMENSION ( : ) :: sum_igr
REAL ( dbl ), ALLOCATABLE, DIMENSION ( : ) :: gauss
REAL ( dbl ), ALLOCATABLE, DIMENSION ( :, : ) :: gdebug
REAL ( dbl ), ALLOCATABLE, DIMENSION ( : ) :: r_delta
INTEGER :: ig, i, idim, natoms, ngtot, isos, lp, mp, np
REAL ( dbl ) :: alpha, epsilon0, ep, em, charge
INTEGER :: gpt, i, idim, natoms, ngtot, isos, lp, mp, np
!------------------------------------------------------------------------------
! allocating
ngtot = pw_grid % ngpts_cut
ALLOCATE (gdebug(3,ngtot),STAT=isos)
IF ( isos /= 0 ) &
CALL stop_memory ( 'force_recip_numer', 'gdebug', 3 * ngtot )
ALLOCATE (sum_igr(ngtot),STAT=isos)
IF ( isos /= 0 ) CALL stop_memory ( 'force_recip_numer', 'sum_igr', ngtot )
ALLOCATE (gauss(ngtot),STAT=isos)
IF ( isos /= 0 ) CALL stop_memory ( 'force_recip_numer', 'gauss', ngtot )
ALLOCATE (r_delta(3),STAT=isos)
IF ( isos /= 0 ) CALL stop_memory ( 'force_recip_numer', 'r_delta', 3 )
! computing the g-vectors
DO ig = 1, ngtot
lp = pw_grid % mapl % pos ( pw_grid % g_hat(1,ig) )
mp = pw_grid % mapm % pos ( pw_grid % g_hat(2,ig) )
np = pw_grid % mapn % pos ( pw_grid % g_hat(3,ig) )
gdebug(1,ig) = 2.0_dbl*pi*(box % h_inv(1,1)*lp+box % h_inv(2,1)* mp &
+box % h_inv(3,1)*np)
gdebug(2,ig) = 2.0_dbl*pi*(box % h_inv(1,2)*lp+box % h_inv(2,2)* mp &
+box % h_inv(3,2)*np)
gdebug(3,ig) = 2.0_dbl*pi*(box % h_inv(1,3)*lp+box % h_inv(2,3)* mp &
+box % h_inv(3,3)*np)
END DO
! defining alpha and epsilon0 (see ewald.f for details)
! defining alpha and epsilon0
alpha = ewald_param % alpha
epsilon0 = ewald_param % eps0
! first initialize the arrays gauss and sum_igr
CALL potential_g_numer ( ep, pnode, sum_igr, gauss, alpha, gdebug, ngtot )
CALL potential_g_numer ( ep, pnode, sum_igr, gauss, alpha, &
pw_grid % g, ngtot )
! initializing numerical force
f_numer = 0.0_dbl
f_numer( :, : ) = 0.0_dbl
! computing the numerical force on each atom
natoms = size(pnode)
@ -478,21 +497,19 @@ SUBROUTINE force_recip_numer ( ewald_param, pnode, box, potparm, &
r_delta ( : ) = pnode(i) % p % r ( : )
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, &
gdebug,ngtot)
CALL de_g_numer(ep, sum_igr, gauss, pnode(i) % p % r , &
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, &
gdebug,ngtot)
f_numer(idim,i) = (em-ep)/epsilon0/box % deth
CALL de_g_numer(em, sum_igr, gauss, pnode(i) % p % r , &
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
END DO
f_numer = f_numer / ( 2.0_dbl * numerical_shift )
f_numer = f_numer / 2.0_dbl / numerical_shift
DEALLOCATE (r_delta,STAT=isos)
IF ( isos /= 0 ) CALL stop_memory ( 'force_recip_numer','r_delta')
DEALLOCATE (gdebug,STAT=isos)
IF ( isos /= 0 ) CALL stop_memory ( 'force_recip_numer','gdebug')
DEALLOCATE (sum_igr,STAT=isos)
IF ( isos /= 0 ) CALL stop_memory ( 'force_recip_numer','sum_igr')
DEALLOCATE (gauss,STAT=isos)
@ -681,13 +698,14 @@ END SUBROUTINE pvbend_numer
!******************************************************************************
SUBROUTINE ptens_numer(pnode,box,potparm,pv_test,delta)
SUBROUTINE ptens_numer(ewald_param,pnode,box,potparm,pv_test,delta)
! computes the numerical pressure tensor
IMPLICIT NONE
! Arguments
TYPE ( particle_node_type ), DIMENSION ( : ), INTENT ( IN ) :: pnode
TYPE ( ewald_parameters_type ), INTENT ( IN ) :: ewald_param
TYPE ( cell_type ), INTENT ( IN ) :: box
REAL ( dbl ), INTENT ( IN ) :: delta
REAL ( dbl ), DIMENSION ( 3, 3 ), INTENT ( OUT ) :: pv_test
@ -731,7 +749,7 @@ SUBROUTINE ptens_numer(pnode,box,potparm,pv_test,delta)
DO i = 1, natoms
x(:,i) = matmul(box_local % hmat,s(:,i))
END DO
CALL getv_ptens(pnode,x,vnb,box_local,potparm)
CALL getv_ptens(ewald_param,pnode,x,vnb,box_local,potparm)
vp = vnb
! tweak again
@ -740,7 +758,7 @@ SUBROUTINE ptens_numer(pnode,box,potparm,pv_test,delta)
DO i = 1, natoms
x(:,i) = matmul(box_local % hmat,s(:,i))
END DO
CALL getv_ptens(pnode,x,vnb,box_local,potparm)
CALL getv_ptens(ewald_param,pnode,x,vnb,box_local,potparm)
vm = vnb
! calculate the derivative
@ -773,7 +791,7 @@ END SUBROUTINE ptens_numer
!******************************************************************************
SUBROUTINE pvg_numer(ewald_param,pnode,box,potparm,pw_grid,pv_test,delta)
SUBROUTINE pvg_numer(ewald_param,pnode,box,pw_grid,pv_test,delta)
! computes the numerical pressure tensor
@ -786,22 +804,21 @@ SUBROUTINE pvg_numer(ewald_param,pnode,box,potparm,pw_grid,pv_test,delta)
TYPE ( cell_type ), INTENT ( IN ) :: box
REAL ( dbl ), INTENT ( IN ) :: delta
REAL ( dbl ), DIMENSION ( :, : ), INTENT ( OUT ) :: pv_test
TYPE ( potentialparm_type ), DIMENSION ( :, : ), INTENT ( IN ) :: potparm
! Locals
TYPE ( pw_grid_type ) :: pw_grid_local
TYPE ( cell_type ) :: box_local
TYPE ( particle_node_type ), DIMENSION ( : ), ALLOCATABLE :: pnode_local
REAL ( dbl ), DIMENSION (3,3) :: dvdh
REAL ( dbl ) :: idelta, vm, vp
REAL ( dbl ), DIMENSION ( :, : ), ALLOCATABLE :: glocal
TYPE ( particle_node_type ), DIMENSION ( : ), ALLOCATABLE :: pnode_local
REAL ( dbl ), DIMENSION ( :, : ), ALLOCATABLE :: s
REAL ( dbl ) :: alpha, epsilon0
COMPLEX ( dbl ), DIMENSION ( : ), ALLOCATABLE :: sum_igr
REAL ( dbl ) :: alpha, epsilon0, cutoff
REAL ( dbl ), DIMENSION ( : ), ALLOCATABLE :: gauss
INTEGER :: i, j, k, ii, jj, ig, ngtot, natoms, isos, lp, mp, np
COMPLEX ( dbl ), DIMENSION ( : ), ALLOCATABLE :: sum_igr
INTEGER :: i, j, k, ii, jj, natoms, isos, gmax, npts_s(3), ngtot
!------------------------------------------------------------------------------
! allocating
natoms = SIZE ( pnode )
ALLOCATE ( s ( 3, natoms ), STAT = isos )
@ -811,14 +828,13 @@ SUBROUTINE pvg_numer(ewald_param,pnode,box,potparm,pw_grid,pv_test,delta)
IF ( isos /= 0 ) CALL stop_memory ( 'pvg_numer', 'sum_igr', ngtot )
ALLOCATE (gauss(ngtot),STAT=isos)
IF ( isos /= 0 ) CALL stop_memory ( 'pvg_numer', 'gauss', ngtot )
ALLOCATE (glocal(3,ngtot),STAT=isos)
IF ( isos /= 0 ) CALL stop_memory ( 'pvg_numer', 'glocal', 3 * ngtot )
ALLOCATE (pnode_local(natoms),STAT=isos)
IF ( isos /= 0 ) CALL stop_memory ( 'pvg_numer', 'pnode_local', natoms )
! assigning the local variables
box_local = box
pnode_local = pnode
pw_grid_local = pw_grid
DO i = 1, natoms
s(:,i) = matmul(box % h_inv,pnode(i) % p % r)
END DO
@ -829,6 +845,10 @@ SUBROUTINE pvg_numer(ewald_param,pnode,box,potparm,pw_grid,pv_test,delta)
! Initializing pv_test
pv_test = 0.0_dbl
gmax = ewald_param % gmax
npts_s = (/ gmax, gmax, gmax /)
CALL pw_find_cutoff ( npts_s, box, cutoff )
! Defining the increments
idelta = 1.0_dbl/(2.0_dbl*delta)
@ -845,27 +865,10 @@ SUBROUTINE pvg_numer(ewald_param,pnode,box,potparm,pw_grid,pv_test,delta)
END DO
! compute the g-vectors
DO ig = 1, ngtot
lp = pw_grid % mapl % pos ( pw_grid % g_hat(1,ig) )
mp = pw_grid % mapm % pos ( pw_grid % g_hat(2,ig) )
np = pw_grid % mapn % pos ( pw_grid % g_hat(3,ig) )
glocal(1,ig) = 2.0_dbl * pi * &
( box % h_inv(1,1) * lp &
+ box % h_inv(2,1) * mp &
+ box % h_inv(3,1) * np )
glocal(2,ig) = 2.0_dbl * pi * &
( box % h_inv(1,2) * lp &
+ box % h_inv(2,2) * mp &
+ box % h_inv(3,2) * np )
glocal(3,ig) = 2.0_dbl * pi * &
( box % h_inv(1,3) * lp &
+ box % h_inv(2,3) * mp &
+ box % h_inv(3,3) * np )
END DO
CALL pw_grid_setup( box_local, pw_grid_local, cutoff)
CALL potential_g_numer ( vp, pnode_local, sum_igr, gauss, alpha, &
glocal, ngtot )
pw_grid_local % g, ngtot )
vp = vp / ( epsilon0 * box_local % deth )
! tweak again
@ -880,26 +883,10 @@ SUBROUTINE pvg_numer(ewald_param,pnode,box,potparm,pw_grid,pv_test,delta)
END DO
! compute the g-vectors
DO ig = 1, ngtot
lp = pw_grid % mapl % pos ( pw_grid % g_hat(1,ig) )
mp = pw_grid % mapm % pos ( pw_grid % g_hat(2,ig) )
np = pw_grid % mapn % pos ( pw_grid % g_hat(3,ig) )
glocal ( 1, ig ) = 2.0_dbl * pi * &
( box % h_inv ( 1, 1 ) * lp &
+ box % h_inv ( 2, 1 ) * mp &
+ box % h_inv ( 3, 1 ) * np )
glocal ( 2, ig ) = 2.0_dbl * pi * &
( box % h_inv ( 1, 2 ) * lp &
+ box % h_inv ( 2, 2 ) * mp &
+ box % h_inv ( 3, 2 ) * np )
glocal ( 3, ig ) = 2.0_dbl * pi * &
( box % h_inv ( 1, 3 ) * lp &
+ box % h_inv ( 2, 3 ) * mp &
+ box % h_inv ( 3, 3 ) * np )
END DO
CALL pw_grid_setup( box_local, pw_grid_local, cutoff)
CALL potential_g_numer ( vm, pnode_local, sum_igr, gauss, alpha, &
glocal, ngtot )
pw_grid_local % g, ngtot )
vm = vm / epsilon0 / box_local % deth
! calculate the derivative
@ -925,10 +912,6 @@ SUBROUTINE pvg_numer(ewald_param,pnode,box,potparm,pw_grid,pv_test,delta)
IF ( isos /= 0 ) CALL stop_memory( 'pvg_numer', 'sum_igr' )
DEALLOCATE (gauss,STAT=isos)
IF ( isos /= 0 ) CALL stop_memory( 'pvg_numer', 'gauss' )
DEALLOCATE (glocal,STAT=isos)
IF ( isos /= 0 ) CALL stop_memory( 'pvg_numer', 'glocal' )
DEALLOCATE (pnode_local,STAT=isos)
IF ( isos /= 0 ) CALL stop_memory( 'pvg_numer', 'pnode_local' )
END SUBROUTINE pvg_numer
@ -959,16 +942,18 @@ SUBROUTINE potential_g_numer(energy,pnode,sum_igr,gauss,alpha, &
energy = 0.0_dbl
sum_igr = 0.0_dbl
DO ig = 1, ngtot
gsq = dot_product(glocal(:,ig),glocal(:,ig))
gauss(ig) = exp(-gsq*.25_dbl/alpha/alpha)/gsq
gsq = DOT_PRODUCT(glocal(:,ig),glocal(:,ig))
IF ( gsq <= 1.0E-10_dbl ) CYCLE
gauss(ig) = exp(-gsq*0.25_dbl/alpha/alpha)/gsq
DO iatom = 1, natoms
gdotr = dot_product(pnode(iatom) % p % r ( : ),glocal(:,ig))
gdotr = DOT_PRODUCT(pnode(iatom) % p % r ( : ), glocal(:,ig))
charge = pnode(iatom) % p % prop % charge
sum_igr(ig) = sum_igr(ig) + charge*cmplx(cos(gdotr),sin(gdotr))
sum_igr(ig) = sum_igr(ig) + charge*CMPLX(COS(gdotr),SIN(gdotr))
END DO
! computing the potential energy
energy = energy + gauss(ig)*sum_igr(ig)*conjg(sum_igr(ig))
energy = energy + gauss(ig)* REAL ( sum_igr ( ig ) * &
CONJG( sum_igr ( ig ) ) )
END DO
END SUBROUTINE potential_g_numer
@ -991,74 +976,60 @@ SUBROUTINE de_g_numer(energy,sum_igr,gauss,r,r_delta,charge,glocal, &
! Locals
INTEGER :: ig
REAL ( dbl ) :: gdotr, gdotr_delta
REAL ( dbl ) :: gdotr, gdotr_delta, gsq
COMPLEX ( dbl ) :: sum
!------------------------------------------------------------------------------
! initialize energy
energy = 0._dbl
sum = ( 0._dbl, 0._dbl )
DO ig = 1, igtot
gsq = DOT_PRODUCT( glocal(:,ig), glocal(:,ig))
IF ( gsq <= 1.0E-10_dbl ) CYCLE
! compute g.r and g.(r+delta)
gdotr = dot_product(r ( : ),glocal(:,ig))
gdotr_delta = dot_product(r_delta ( : ),glocal(:,ig))
gdotr = DOT_PRODUCT ( r ( : ), glocal ( :, ig ) )
gdotr_delta = DOT_PRODUCT ( r_delta ( : ), glocal ( :, ig ) )
! subtract off exp(ig.r) and add exp(ig.(r+delta))
sum = sum_igr(ig) - charge*cmplx(cos(gdotr),sin(gdotr)) + &
charge*cmplx(cos(gdotr_delta),sin(gdotr_delta))
sum = sum_igr ( ig ) - charge * CMPLX ( COS( gdotr ), SIN( gdotr ) ) + &
charge * CMPLX( COS( gdotr_delta ), SIN( gdotr_delta ) )
! recompute energy
energy = energy + gauss(ig)*sum*conjg(sum)
energy = energy + gauss(ig) * REAL ( sum * CONJG ( sum ) )
END DO
END SUBROUTINE de_g_numer
!******************************************************************************
SUBROUTINE energy_recip_numer ( ewald_param, pnode, box, potparm, pw_grid, &
SUBROUTINE energy_recip_numer ( ewald_param, pnode, box, pw_grid, &
energy_numer, ngtot )
IMPLICIT NONE
! Arguments
TYPE ( ewald_parameters_type ), INTENT ( IN ) :: ewald_param
TYPE ( pw_grid_type ), INTENT ( IN ) :: pw_grid
TYPE ( pw_grid_type ), INTENT ( IN ):: pw_grid
TYPE ( particle_node_type ), DIMENSION ( : ), INTENT ( IN ) :: pnode
TYPE ( cell_type ), INTENT ( IN ) :: box
TYPE ( potentialparm_type ), DIMENSION ( :, : ), INTENT ( IN ) :: potparm
REAL ( dbl ), INTENT ( OUT ) :: energy_numer
INTEGER, INTENT ( IN ) :: ngtot
! Locals
COMPLEX ( dbl ), ALLOCATABLE, DIMENSION ( : ) :: sum_igr
REAL ( dbl ), ALLOCATABLE, DIMENSION ( : ) :: gauss
REAL ( dbl ), ALLOCATABLE, DIMENSION ( :, : ) :: gdebug
REAL ( dbl ) :: alpha, epsilon0
INTEGER :: ig, i, idim, isos, lp, mp, np
REAL ( dbl ) :: alpha, epsilon0, ep, em, charge
!------------------------------------------------------------------------------
! allocating
ALLOCATE (gdebug(3,ngtot),STAT=isos)
IF ( isos /= 0 ) &
CALL stop_memory ( 'energy_recip_numer', 'gdebug', 3 * ngtot )
ALLOCATE ( sum_igr ( ngtot ), STAT = isos )
IF ( isos /= 0 ) CALL stop_memory ( 'energy_recip_numer', 'sum_igr', ngtot )
ALLOCATE (gauss(ngtot),STAT=isos)
ALLOCATE ( gauss ( ngtot ), STAT = isos )
IF ( isos /= 0 ) CALL stop_memory ( 'energy_recip_numer', 'gauss', ngtot )
! computing the g-vectors
DO ig = 1, ngtot
lp = pw_grid % mapl % pos ( pw_grid % g_hat(1,ig) )
mp = pw_grid % mapm % pos ( pw_grid % g_hat(2,ig) )
np = pw_grid % mapn % pos ( pw_grid % g_hat(3,ig) )
gdebug(1,ig) = 2.0_dbl*pi*(box % h_inv(1,1)*lp+box % h_inv(2,1)*mp &
+box % h_inv(3,1)*np)
gdebug(2,ig) = 2.0_dbl*pi*(box % h_inv(1,2)*lp+box % h_inv(2,2)*mp &
+box % h_inv(3,2)*np)
gdebug(3,ig) = 2.0_dbl*pi*(box % h_inv(1,3)*lp+box % h_inv(2,3)*mp &
+box % h_inv(3,3)*np)
END DO
! defining alpha and epsilon0 (see ewald.f for details)
! defining alpha and epsilon0
alpha = ewald_param % alpha
epsilon0 = ewald_param % eps0
@ -1067,11 +1038,8 @@ SUBROUTINE energy_recip_numer ( ewald_param, pnode, box, potparm, pw_grid, &
! computing the energy
CALL potential_g_numer(energy_numer,pnode,sum_igr,gauss,alpha, &
gdebug,ngtot)
pw_grid % g, ngtot)
energy_numer = energy_numer/epsilon0/box % deth
DEALLOCATE (gdebug,STAT=isos)
IF ( isos /= 0 ) CALL stop_memory ( 'energy_recip_numer', 'gdebug' )
DEALLOCATE (sum_igr,STAT=isos)
IF ( isos /= 0 ) CALL stop_memory ( 'energy_recip_numer', 'sum_igr' )
DEALLOCATE (gauss,STAT=isos)
@ -1157,7 +1125,7 @@ END SUBROUTINE vnumer_bends
!******************************************************************************
SUBROUTINE getv_ptens ( pnode, x, vnb, box, potparm )
SUBROUTINE getv_ptens ( ewald_param, pnode, x, vnb, box, potparm )
IMPLICIT NONE
@ -1167,6 +1135,7 @@ SUBROUTINE getv_ptens ( pnode, x, vnb, box, potparm )
TYPE (cell_type), INTENT ( INOUT ) :: box
TYPE (particle_node_type ), INTENT ( IN ), DIMENSION ( : ) :: pnode
TYPE (potentialparm_type), INTENT ( IN ), DIMENSION ( :, : ) :: potparm
TYPE (ewald_parameters_type), INTENT (IN) :: ewald_param
! Locals
INTEGER :: i, j, ii, jj, iatomtype, jatomtype, iexclude, natoms, id
@ -1248,9 +1217,13 @@ SUBROUTINE getv_ptens ( pnode, x, vnb, box, potparm )
END DO EXCL
NOTMATCH: IF ( .NOT. match) THEN
CALL potential_f ( rijsq, potparm, qi, qj, &
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
!
! summing up the potential energy
!
@ -1279,9 +1252,14 @@ SUBROUTINE getv_ptens ( pnode, x, vnb, box, potparm )
CALL find_image(s,perd,vec,box % hmat,rijsq,rij)
IF ( rijsq <= potparm ( iatomtype, jatomtype ) % rcutsq ) &
THEN
CALL potential_f ( rijsq, potparm, qi, qj, &
iatomtype,jatomtype, energy )
vnb = vnb + energy
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
vnb = vnb + energy
END IF
END DO
END DO
@ -1302,9 +1280,14 @@ SUBROUTINE getv_ptens ( pnode, x, vnb, box, potparm )
CALL find_image ( s, perd, vec, box % hmat, rijsq, rij )
IF ( rijsq <= potparm ( iatomtype, iatomtype ) % rcutsq ) &
THEN
CALL potential_f ( rijsq, potparm, qi, qi, &
iatomtype, iatomtype, energy )
vnb = vnb + energy
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
vnb = vnb + energy
END IF
END DO

7
src/fist_intra_force.F
View file

@ -10,6 +10,7 @@ MODULE fist_intra_force
get_pv_bond, get_pv_bend, get_pv_torsion
USE molecule_types, ONLY : molecule_structure_type, linklist_bonds, &
linklist_bends
USE timings, ONLY : timeset, timestop
IMPLICIT NONE
@ -39,9 +40,11 @@ SUBROUTINE force_intra_control ( molecule, v_bond, v_bend, pvbond, pvbend, &
REAL ( dbl ), DIMENSION (3) :: rij, b12, b32, g1, g2, g3
TYPE (linklist_bonds), POINTER :: llbond
TYPE (linklist_bends), POINTER :: llbend
INTEGER :: ibond, ibend, imol
INTEGER :: ibond, ibend, imol, handle
!------------------------------------------------------------------------------
CALL timeset ( 'FORCE_INTRA_CONTROL','I',' ',handle )
IF ( PRESENT ( f_bond)) f_bond = 0._dbl
IF ( PRESENT ( f_bend)) f_bend = 0._dbl
@ -110,6 +113,8 @@ SUBROUTINE force_intra_control ( molecule, v_bond, v_bend, pvbond, pvbend, &
END DO BEND
END DO MOL
CALL timestop ( 0._dbl, handle )
END SUBROUTINE force_intra_control

19
src/force_control.F
View file

@ -7,14 +7,12 @@ MODULE force_control
USE ewald_parameters_types, ONLY : ewald_parameters_type
USE fist_force, ONLY : fist_force_control => force_control
USE fist_global, ONLY : fistpar
USE global_types, ONLY : global_environment_type
USE tbmd_force, ONLY : tbmd_force_control => force_control
USE tbmd_global, ONLY : tbmdpar
USE kinds, ONLY : dbl
USE md, ONLY : simulation_parameters_type, thermodynamic_type
USE pair_potential, ONLY : potentialparm_type
USE stop_program, ONLY : stop_prg
USE structure_types, ONLY : structure_type
USE structure_types, ONLY : structure_type, interaction_type
IMPLICIT NONE
@ -25,7 +23,7 @@ CONTAINS
!******************************************************************************
SUBROUTINE force ( struc, potparm, thermo, simpar, ewald_param )
SUBROUTINE force ( struc, inter, thermo, simpar, ewald_param, globenv )
IMPLICIT NONE
@ -34,7 +32,8 @@ SUBROUTINE force ( struc, potparm, thermo, simpar, ewald_param )
TYPE ( thermodynamic_type ), INTENT ( INOUT ) :: thermo
TYPE ( simulation_parameters_type ), INTENT ( IN ) :: simpar
TYPE ( ewald_parameters_type ), INTENT ( INOUT ) :: ewald_param
TYPE ( potentialparm_type ), DIMENSION ( :, : ), INTENT ( IN ) :: potparm
TYPE ( interaction_type ), INTENT ( IN ) :: inter
TYPE ( global_environment_type), INTENT ( IN ) :: globenv
!------------------------------------------------------------------------------
@ -44,13 +43,13 @@ SUBROUTINE force ( struc, potparm, thermo, simpar, ewald_param )
CASE ( 'FIST' )
CALL fist_force_control ( struc % molecule, struc % pnode, struc % part, &
struc % box, thermo, potparm, ewald_param, simpar % ensemble, &
fistpar )
struc % box, thermo, inter%potparm, ewald_param, simpar % ensemble, &
globenv )
CASE ( 'TBMD' )
CALL tbmd_force_control ( struc % molecule, struc % pnode, struc % part, &
struc % box, thermo, potparm, ewald_param, simpar % ensemble, &
tbmdpar )
struc % box, thermo, inter%potparm, ewald_param, simpar % ensemble, &
globenv )
END SELECT

91
src/integrator.F
View file

@ -9,6 +9,7 @@ MODULE integrator
! pv_constraint, initialize_roll, rattle_roll_control, &
! shake_roll_control
USE constraint
USE global_types, ONLY: global_environment_type
USE eigenvalueproblems, ONLY : diagonalise
USE ewald_parameters_types, ONLY : ewald_parameters_type
USE force_control, ONLY : force
@ -22,11 +23,10 @@ MODULE integrator
linklist_atoms
USE message_passing, ONLY : mp_sum
USE nose, ONLY : extended_parameters_type, lnhc, lnhcp, lnhcpf
USE pair_potential, ONLY : potentialparm_type
USE particle_types, ONLY : particle_type
USE simulation_cell, ONLY : get_cell_param
USE stop_program, ONLY : stop_prg, stop_memory
USE structure_types, ONLY : structure_type
USE structure_types, ONLY : structure_type, interaction_type
USE timings, ONLY : timeset, timestop
USE util, ONLY : get_unit
@ -65,11 +65,15 @@ MODULE integrator
INTEGER :: crd, vel, ptn, ene, tem, scr
INTEGER :: icrd, ivel, iptens, iener, itemp, idump, iscreen
TYPE ( global_environment_type ) :: intenv
!******************************************************************************
CONTAINS
!******************************************************************************
SUBROUTINE velocity_verlet ( itimes, constant, simpar, potparm, &
SUBROUTINE velocity_verlet ( itimes, constant, simpar, inter, &
thermo, struc, ewald_param, nhcp )
IMPLICIT NONE
@ -82,26 +86,30 @@ SUBROUTINE velocity_verlet ( itimes, constant, simpar, potparm, &
TYPE ( structure_type ), INTENT ( INOUT ) :: struc
TYPE ( extended_parameters_type ), INTENT ( INOUT ) :: nhcp
TYPE ( ewald_parameters_type ), INTENT ( INOUT ) :: ewald_param
TYPE ( potentialparm_type ), DIMENSION ( :, : ), INTENT ( IN ) :: potparm
TYPE ( interaction_type ), INTENT ( IN ) :: inter
INTEGER :: handle
!------------------------------------------------------------------------------
CALL timeset ( 'VERLET', 'I', ' ', handle )
SELECT CASE ( simpar % ensemble )
CASE DEFAULT
CALL stop_prg ( 'velocity_verlet','integrator not implemented')
CASE ( 'NVE')
CALL nve(itimes,constant,simpar,potparm,thermo,struc,ewald_param)
CALL nve(itimes,constant,simpar,inter,thermo,struc,ewald_param)
CALL energy(itimes,constant,simpar,struc,thermo,nhcp)
CASE ( 'NVT')
CALL nvt(itimes,constant,simpar,potparm,thermo,struc,ewald_param,nhcp)
CALL nvt(itimes,constant,simpar,inter,thermo,struc,ewald_param,nhcp)
CALL energy(itimes,constant,simpar,struc,thermo,nhcp)
CASE ( 'NPT_I')
CALL npt_i(itimes,constant,simpar,potparm,thermo,struc,ewald_param,nhcp)
CALL npt_i(itimes,constant,simpar,inter,thermo,struc,ewald_param,nhcp)
CALL energy(itimes,constant,simpar,struc,thermo,nhcp)
CASE ( 'NPT_F')
CALL npt_f(itimes,constant,simpar,potparm,thermo,struc,ewald_param,nhcp)
CALL npt_f(itimes,constant,simpar,inter,thermo,struc,ewald_param,nhcp)
CALL energy(itimes,constant,simpar,struc,thermo,nhcp)
END SELECT
CALL timestop ( zero, handle )
END SUBROUTINE velocity_verlet
@ -110,20 +118,20 @@ END SUBROUTINE velocity_verlet
! 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 ( ion, group, iw, mdio )
SUBROUTINE set_integrator ( globenv, mdio )
IMPLICIT NONE
! Arguments
LOGICAL, INTENT (IN) :: ion
INTEGER, INTENT (IN) :: group
INTEGER, INTENT (IN) :: iw
TYPE(mdio_parameters_type), INTENT (IN) :: mdio
TYPE ( global_environment_type ), INTENT (IN) :: globenv
!------------------------------------------------------------------------------
ionode = ion
int_group = group
scr = iw
intenv = globenv
ionode = globenv % ionode
int_group = globenv % group
scr = globenv % scr
crd_file_name = mdio % crd_file_name
vel_file_name = mdio % vel_file_name
@ -140,12 +148,6 @@ SUBROUTINE set_integrator ( ion, group, iw, mdio )
idump = mdio % idump
iscreen = mdio % iscreen
crd = get_unit()
vel = get_unit()
ptn = get_unit()
ene = get_unit()
tem = get_unit()
END SUBROUTINE set_integrator
!******************************************************************************
@ -197,13 +199,18 @@ SUBROUTINE energy ( itimes, constant, simpar, struc, thermo, nhcp )
!------------------------------------------------------------------------------
CALL timeset ( 'ENERGY', 'I', ' ', handle )
CALL timeset ( 'ENERGY', 'E', ' ', handle )
IF ( ionode .AND. itimes == 0 ) THEN
tem = get_unit()
OPEN ( UNIT = tem, FILE = temp_file_name )
ene = get_unit()
OPEN ( UNIT = ene, FILE = ener_file_name )
crd = get_unit()
OPEN ( UNIT = crd, FILE = crd_file_name )
vel = get_unit()
OPEN ( UNIT = vel, FILE = vel_file_name )
ptn = get_unit()
OPEN ( UNIT = ptn, FILE = ptens_file_name )
END IF
@ -444,7 +451,7 @@ END SUBROUTINE energy
! nve integrator for particle positions & momenta
! velocity Verlet
SUBROUTINE nve(itimes,constant,simpar,potparm,thermo,struc,ewald_param)
SUBROUTINE nve(itimes,constant,simpar,inter,thermo,struc,ewald_param)
IMPLICIT NONE
INTEGER, INTENT ( IN ) :: itimes
REAL ( dbl ), INTENT ( INOUT ) :: constant
@ -452,7 +459,7 @@ SUBROUTINE nve(itimes,constant,simpar,potparm,thermo,struc,ewald_param)
TYPE (thermodynamic_type ), INTENT ( INOUT ) :: thermo
TYPE (structure_type ), INTENT ( INOUT ) :: struc
TYPE (ewald_parameters_type ), INTENT ( INOUT ) :: ewald_param
TYPE (potentialparm_type ), INTENT ( IN ), DIMENSION ( :, : ) :: potparm
TYPE (interaction_type ), INTENT ( IN ) :: inter
INTEGER :: i, nnodes
REAL ( dbl ) :: dtom
@ -489,7 +496,7 @@ SUBROUTINE nve(itimes,constant,simpar,potparm,thermo,struc,ewald_param)
!
! get new forces
!
CALL force ( struc, potparm, thermo, simpar, ewald_param )
CALL force ( struc, inter, thermo, simpar, ewald_param, intenv )
!
! second half of velocity verlet
@ -519,7 +526,7 @@ END SUBROUTINE nve
!******************************************************************************
SUBROUTINE nvt(itimes,constant,simpar,potparm,thermo,struc, &
SUBROUTINE nvt(itimes,constant,simpar,inter,thermo,struc, &
ewald_param,nhcp)
IMPLICIT NONE
@ -532,7 +539,7 @@ SUBROUTINE nvt(itimes,constant,simpar,potparm,thermo,struc, &
TYPE (structure_type ), INTENT ( INOUT ) :: struc
TYPE (extended_parameters_type ), INTENT ( INOUT ) :: nhcp
TYPE (ewald_parameters_type ), INTENT ( INOUT ) :: ewald_param
TYPE (potentialparm_type ), INTENT ( IN ), DIMENSION ( :, : ) :: potparm
TYPE (interaction_type ), INTENT ( IN ) :: inter
! Locals
INTEGER :: i, nnodes
@ -572,7 +579,7 @@ SUBROUTINE nvt(itimes,constant,simpar,potparm,thermo,struc, &
!
! get new forces
!
CALL force ( struc, potparm, thermo, simpar, ewald_param )
CALL force ( struc, inter, thermo, simpar, ewald_param, intenv )
!
! second half of velocity verlet
@ -604,7 +611,7 @@ END SUBROUTINE nvt
!******************************************************************************
SUBROUTINE npt_i ( itimes, constant, simpar, potparm, thermo, struc, &
SUBROUTINE npt_i ( itimes, constant, simpar, inter, thermo, struc, &
ewald_param, nhcp )
IMPLICIT NONE
@ -617,7 +624,7 @@ SUBROUTINE npt_i ( itimes, constant, simpar, potparm, thermo, struc, &
TYPE (structure_type ), INTENT ( INOUT ) :: struc
TYPE (extended_parameters_type ), INTENT ( INOUT ) :: nhcp
TYPE (ewald_parameters_type ), INTENT ( INOUT ) :: ewald_param
TYPE (potentialparm_type ), INTENT ( IN ), DIMENSION ( :, : ) :: potparm
TYPE (interaction_type ), INTENT ( IN ) :: inter
! Locals
INTEGER :: i, nnodes, iroll
@ -700,7 +707,7 @@ SUBROUTINE npt_i ( itimes, constant, simpar, potparm, thermo, struc, &
!
! get new forces
!
CALL force(struc,potparm,thermo,simpar,ewald_param)
CALL force(struc,inter,thermo,simpar,ewald_param,intenv)
!
! second half of velocity verlet
@ -757,10 +764,13 @@ SUBROUTINE pressure ( pnode, thermo )
TYPE ( thermodynamic_type ), INTENT ( INOUT ) :: thermo
! Locals
INTEGER :: i, j, iatom, nnodes
INTEGER :: i, j, iatom, nnodes, handle
REAL ( dbl ) :: mfl
!------------------------------------------------------------------------------
CALL timeset ( 'PRESSURE', 'E', 'Mflops', handle )
thermo%pv_kin = 0._dbl
nnodes = size(pnode)
DO i = 1, 3
@ -772,6 +782,7 @@ SUBROUTINE pressure ( pnode, thermo )
END DO
END DO
END DO
mfl = REAL( 9 * nnodes, dbl ) * 2._dbl * 1.e-6_dbl
#if defined(__parallel)
CALL mp_sum(thermo%pv_kin,int_group)
@ -780,6 +791,8 @@ SUBROUTINE pressure ( pnode, thermo )
! total virial
thermo%ptens = thermo%pv + thermo%pv_kin + thermo%pv_const
CALL timestop ( mfl, handle )
END SUBROUTINE pressure
!******************************************************************************
@ -794,11 +807,13 @@ SUBROUTINE update_structure(struc,task)
! Locals
REAL ( dbl ), ALLOCATABLE, DIMENSION ( :, : ) :: atot
INTEGER :: natoms, ios, imol, atombase, iat, i
INTEGER :: natoms, ios, imol, atombase, iat, i, handle
TYPE (linklist_atoms), POINTER :: current_atom
!------------------------------------------------------------------------------
CALL timeset ( 'UPDATE_STRUCTURE', 'E', ' ', handle )
natoms = size(struc%part)
ALLOCATE (atot(3,natoms),STAT=ios)
IF ( ios /= 0 ) CALL stop_memory ( 'update_structure', 'atot', 3 * natoms )
@ -840,6 +855,8 @@ SUBROUTINE update_structure(struc,task)
DEALLOCATE ( atot, STAT = ios )
IF ( ios /= 0 ) CALL stop_memory ( 'update_structure', 'atot' )
CALL timestop ( zero, handle )
END SUBROUTINE update_structure
@ -910,7 +927,7 @@ END SUBROUTINE set
!******************************************************************************
SUBROUTINE npt_f ( itimes, constant, simpar, potparm, thermo, struc, &
SUBROUTINE npt_f ( itimes, constant, simpar, inter, thermo, struc, &
ewald_param, nhcp )
IMPLICIT NONE
@ -923,7 +940,7 @@ SUBROUTINE npt_f ( itimes, constant, simpar, potparm, thermo, struc, &
TYPE (structure_type ), INTENT ( INOUT ) :: struc
TYPE (extended_parameters_type ), INTENT ( INOUT ) :: nhcp
TYPE (ewald_parameters_type ), INTENT ( INOUT ) :: ewald_param
TYPE (potentialparm_type ), INTENT ( IN ), DIMENSION ( :, : ) :: potparm
TYPE (interaction_type ), INTENT ( IN ) :: inter
! Locals
INTEGER :: i, j, nnodes, iroll
@ -1009,7 +1026,7 @@ SUBROUTINE npt_f ( itimes, constant, simpar, potparm, thermo, struc, &
!
! get new forces
!
CALL force ( struc, potparm, thermo, simpar, ewald_param )
CALL force ( struc, inter, thermo, simpar, ewald_param, intenv )
!
! second half of velocity verlet

2773
src/k290.F Normal file
View file

File diff suppressed because it is too large Load diff

64
src/kpoint_initialization.F Normal file
View file

@ -0,0 +1,64 @@
!------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!------------------------------------------------------------------------------!
!
MODULE kpoint_initialization
!
!------------------------------------------------------------------------------!
USE kinds, ONLY : dbl
USE brillouin, ONLY : kpoint_type
USE stop_program, ONLY : stop_memory
USE particle_types, ONLY : particle_type
USE global_types, ONLY : global_environment_type
USE k290, ONLY : set_k290, set_k290_atoms, kp_sym_gen
!
IMPLICIT NONE
!
PRIVATE
!
PUBLIC :: initialize_kpoints
!------------------------------------------------------------------------------!
!
CONTAINS
!
!------------------------------------------------------------------------------!
SUBROUTINE initialize_kpoints(inpar,kp,symm,hmat,part)
IMPLICIT NONE
TYPE (global_environment_type) :: inpar
TYPE (kpoint_type) :: kp
LOGICAL, INTENT (IN) :: symm
REAL (dbl), INTENT (IN) :: hmat(:,:)
TYPE (particle_type), INTENT (IN) :: part(:)
REAL (dbl), ALLOCATABLE :: coor(:,:)
INTEGER, ALLOCATABLE :: types(:)
INTEGER :: i, nat, isos
IF (symm .OR. kp%scheme=='MONKHORST-PACK' .OR. kp%scheme=='MACDONALD') &
THEN
CALL set_k290(hmat,kp%nk,shift=kp%shift,symm=kp%symmetry, &
unit=inpar%scr)
nat = size(part)
ALLOCATE (coor(3,nat),STAT=isos)
IF (isos/=0) CALL stop_memory('initialize_kpoints','coor',3*nat)
ALLOCATE (types(nat),STAT=isos)
IF (isos/=0) CALL stop_memory('initialize_kpoints','types',nat)
DO i = 1, nat
coor(:,i) = part(i) %r
types(i) = part(i) %prop%ptype
END DO
CALL set_k290_atoms(coor,types)
DEALLOCATE (coor,STAT=isos)
IF (isos/=0) CALL stop_memory('initialize_kpoints','coor')
DEALLOCATE (types,STAT=isos)
IF (isos/=0) CALL stop_memory('initialize_kpoints','types')
CALL kp_sym_gen
END IF
END SUBROUTINE initialize_kpoints
!------------------------------------------------------------------------------!
END MODULE kpoint_initialization
!------------------------------------------------------------------------------!

15
src/kpoints.F Normal file
View file

@ -0,0 +1,15 @@
MODULE kpoints
USE kinds, ONLY : dbl
IMPLICIT NONE
TYPE kpts_type
REAL ( dbl ), DIMENSION ( :, : ), POINTER :: xkpt
REAL ( dbl ), DIMENSION ( : ), POINTER :: wkpt
END TYPE kpts_type
TYPE ( kpts_type ) :: kset
END MODULE kpoints

25
src/lib/Makefile Normal file
View file

@ -0,0 +1,25 @@
.SUFFIXES: .o .F
include ../MACHINEDEFS
OBJECTS_fftsg = mltfftsg.o ctrig.o fftstp.o fftpre.o fftrot.o
#################################
LIB1 = fftsg
all: $(LIB1)
$(LIB1): $(OBJECTS_fftsg)
$(AR) libfftsg.a $(OBJECTS_fftsg)
%.o : %.F
$(FC) -c $(FFLAGS) $*.F
clean:
rm -f *.o *.mod
realclean:
rm -f *.a *.o *.mod *~ F*.f *.lst
####################

111
src/lib/ctrig.F Normal file
View file

@ -0,0 +1,111 @@
!-----------------------------------------------------------------------------!
! Copyright by Stefan Goedecker, Lausanne, Switzerland, August 1, 1991
! modified by Stefan Goedecker, Cornell, Ithaca, USA, March 25, 1994
! modified by Stefan Goedecker, Stuttgart, Germany, October 6, 1995
! Commercial use is prohibited
! without the explicit permission of the author.
!-----------------------------------------------------------------------------!
SUBROUTINE ctrig ( n, trig, after, before, now, isign, ic )
IMPLICIT NONE
INTEGER, PARAMETER :: dbl = SELECTED_REAL_KIND ( 14, 200 )
INTEGER, INTENT ( IN ) :: n
INTEGER, INTENT ( IN ) :: isign
INTEGER, INTENT ( OUT ) :: ic
INTEGER, DIMENSION ( 7 ), INTENT ( OUT ) :: after, before, now
REAL ( dbl ) , DIMENSION ( 2, 1024 ), INTENT ( OUT ) :: trig
INTEGER :: i, j, itt
REAL ( dbl ) :: twopi, angle
INTEGER, PARAMETER :: nt = 82
INTEGER, DIMENSION ( 7, nt ), PARAMETER :: idata = RESHAPE ((/ &
3, 3, 1, 1, 1, 1, 1, 4, 4, 1, 1, 1, 1, 1, &
5, 5, 1, 1, 1, 1, 1, 6, 6, 1, 1, 1, 1, 1, &
8, 8, 1, 1, 1, 1, 1, 9, 3, 3, 1, 1, 1, 1, &
12, 4, 3, 1, 1, 1, 1, 15, 5, 3, 1, 1, 1, 1, &
16, 4, 4, 1, 1, 1, 1, 18, 6, 3, 1, 1, 1, 1, &
20, 5, 4, 1, 1, 1, 1, 24, 8, 3, 1, 1, 1, 1, &
25, 5, 5, 1, 1, 1, 1, 27, 3, 3, 3, 1, 1, 1, &
30, 6, 5, 1, 1, 1, 1, 32, 8, 4, 1, 1, 1, 1, &
36, 4, 3, 3, 1, 1, 1, 40, 8, 5, 1, 1, 1, 1, &
45, 5, 3, 3, 1, 1, 1, 48, 4, 4, 3, 1, 1, 1, &
54, 6, 3, 3, 1, 1, 1, 60, 5, 4, 3, 1, 1, 1, &
64, 4, 4, 4, 1, 1, 1, 72, 8, 3, 3, 1, 1, 1, &
75, 5, 5, 3, 1, 1, 1, 80, 5, 4, 4, 1, 1, 1, &
81, 3, 3, 3, 3, 1, 1, 90, 6, 5, 3, 1, 1, 1, &
96, 8, 4, 3, 1, 1, 1, 100, 5, 5, 4, 1, 1, 1, &
108, 4, 3, 3, 3, 1, 1, 120, 8, 5, 3, 1, 1, 1, &
125, 5, 5, 5, 1, 1, 1, 128, 8, 4, 4, 1, 1, 1, &
135, 5, 3, 3, 3, 1, 1, 144, 4, 4, 3, 3, 1, 1, &
150, 6, 5, 5, 1, 1, 1, 160, 8, 5, 4, 1, 1, 1, &
162, 6, 3, 3, 3, 1, 1, 180, 5, 4, 3, 3, 1, 1, &
192, 4, 4, 4, 3, 1, 1, 200, 8, 5, 5, 1, 1, 1, &
216, 8, 3, 3, 3, 1, 1, 225, 5, 5, 3, 3, 1, 1, &
240, 5, 4, 4, 3, 1, 1, 243, 3, 3, 3, 3, 3, 1, &
256, 4, 4, 4, 4, 1, 1, 270, 6, 5, 3, 3, 1, 1, &
288, 8, 4, 3, 3, 1, 1, 300, 5, 5, 4, 3, 1, 1, &
320, 5, 4, 4, 4, 1, 1, 324, 4, 3, 3, 3, 3, 1, &
360, 8, 5, 3, 3, 1, 1, 375, 5, 5, 5, 3, 1, 1, &
384, 8, 4, 4, 3, 1, 1, 400, 5, 5, 4, 4, 1, 1, &
405, 5, 3, 3, 3, 3, 1, 432, 4, 4, 3, 3, 3, 1, &
450, 6, 5, 5, 3, 1, 1, 480, 8, 5, 4, 3, 1, 1, &
486, 6, 3, 3, 3, 3, 1, 500, 5, 5, 5, 4, 1, 1, &
512, 8, 4, 4, 4, 1, 1, 540, 5, 4, 3, 3, 3, 1, &
576, 4, 4, 4, 3, 3, 1, 600, 8, 5, 5, 3, 1, 1, &
625, 5, 5, 5, 5, 1, 1, 640, 8, 5, 4, 4, 1, 1, &
648, 8, 3, 3, 3, 3, 1, 675, 5, 5, 3, 3, 3, 1, &
720, 5, 4, 4, 3, 3, 1, 729, 3, 3, 3, 3, 3, 3, &
750, 6, 5, 5, 5, 1, 1, 768, 4, 4, 4, 4, 3, 1, &
800, 8, 5, 5, 4, 1, 1, 810, 6, 5, 3, 3, 3, 1, &
864, 8, 4, 3, 3, 3, 1, 900, 5, 5, 4, 3, 3, 1, &
960, 5, 4, 4, 4, 3, 1, 972, 4, 3, 3, 3, 3, 3, &
1000, 8, 5, 5, 5, 1, 1, 1024, 4, 4, 4, 4, 4, 1 /),(/7,nt/))
!-----------------------------------------------------------------------------!
mloop: DO i = 1, nt
IF ( n == idata ( 1, i ) ) THEN
ic=0
DO j = 1, 6
itt = idata ( 1 + j, i )
IF ( itt > 1 ) THEN
ic = ic + 1
now ( j ) = idata ( 1 + j, i )
ELSE
EXIT mloop
END IF
END DO
EXIT mloop
END IF
IF ( i == nt ) THEN
WRITE ( *, '(A,i5,A)' ) " Value of ",n, &
" not allowed for fft, allowed values are:"
WRITE ( *, '(15i5)' ) ( idata ( 1, j ), j = 1, nt )
stop 'ctrig'
END IF
END DO mloop
after ( 1 ) = 1
before ( ic ) = 1
DO i = 2, ic
after ( i ) = after ( i - 1 ) * now ( i - 1 )
before ( ic - i + 1 ) = before ( ic - i + 2 ) * now ( ic - i + 2 )
END DO
twopi = 8._dbl * ATAN ( 1._dbl )
angle = isign * twopi / REAL ( n, dbl )
trig ( 1, 1 ) = 1._dbl
trig ( 2, 1 ) = 0._dbl
DO i = 1, n - 1
trig ( 1, i + 1 ) = cos ( REAL ( i, dbl ) * angle )
trig ( 2, i + 1 ) = sin ( REAL ( i, dbl ) * angle )
END DO
!-----------------------------------------------------------------------------!
END SUBROUTINE ctrig
!-----------------------------------------------------------------------------!

1051
src/lib/fftpre.F Normal file
View file

File diff suppressed because it is too large Load diff

1051
src/lib/fftrot.F Normal file
View file

File diff suppressed because it is too large Load diff

1051
src/lib/fftstp.F Normal file
View file

File diff suppressed because it is too large Load diff

133
src/lib/mltfftsg.F Normal file
View file

@ -0,0 +1,133 @@
!-----------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
SUBROUTINE mltfftsg ( transa, transb, a, ldax, lday, b, ldbx, ldby, &
n, m, isign, scale )
IMPLICIT NONE
INTEGER, PARAMETER :: dbl = SELECTED_REAL_KIND ( 14, 200 )
! Arguments
CHARACTER ( LEN = 1 ), INTENT ( IN ) :: transa*1, transb*1
INTEGER, INTENT ( IN ) :: ldax, lday, ldbx, ldby, n, m, isign
COMPLEX ( dbl ), INTENT ( INOUT ) :: a ( ldax, lday ), b ( ldbx, ldby )
REAL ( dbl ), INTENT ( IN ) :: scale
! Variables
INTEGER, SAVE :: ncache
INTEGER :: after ( 20 ), now ( 20 ), before ( 20 )
REAL ( dbl ) :: trig ( 2, 1024 )
LOGICAL :: tscal
COMPLEX ( dbl ), DIMENSION ( :, : ), ALLOCATABLE, SAVE :: z
INTEGER :: length, isig, ic, lot, itr, nfft, i, inzee, j
!-----------------------------------------------------------------------------!
IF( .NOT. ALLOCATED ( Z ) ) THEN
ncache = get_cache_size ( )
LENGTH = 2 * ( NCACHE / 4 + 1 )
ALLOCATE ( Z ( LENGTH, 2 ) )
END IF
ISIG = -ISIGN
TSCAL = ( ABS ( SCALE -1._dbl ) > 1.e-12_dbl )
CALL CTRIG ( N, TRIG, AFTER, BEFORE, NOW, ISIG, IC )
LOT = NCACHE / ( 4 * N )
LOT = LOT - MOD ( LOT + 1, 2 )
LOT = MAX ( 1, LOT )
DO ITR = 1, M, LOT
NFFT = MIN ( M - ITR + 1, LOT )
IF ( TRANSA == 'N' .OR. TRANSA == 'n' ) THEN
CALL FFTPRE ( NFFT, NFFT, LDAX, LOT, N, A ( 1, ITR ), Z ( 1, 1 ), &
TRIG, NOW ( 1 ), AFTER ( 1 ), BEFORE ( 1 ), ISIG )
ELSE
CALL FFTSTP ( LDAX, NFFT, N, LOT, N, A ( ITR, 1 ), Z ( 1, 1 ), &
TRIG, NOW ( 1 ), AFTER ( 1 ), BEFORE ( 1 ), ISIG )
ENDIF
IF ( TSCAL ) THEN
IF ( LOT == NFFT ) THEN
CALL DSCAL ( 2 * LOT * N, SCALE, Z ( 1, 1 ), 1 )
ELSE
DO I = 1, N
CALL DSCAL ( 2 * NFFT, SCALE, Z ( LOT * ( I - 1 ) + 1, 1 ), 1 )
END DO
END IF
END IF
IF(IC.EQ.1) THEN
IF(TRANSB == 'N'.OR.TRANSB == 'n') THEN
CALL ZGETMO(Z(1,1),LOT,NFFT,N,B(1,ITR),LDBX)
ELSE
CALL MATMOV(NFFT,N,Z(1,1),LOT,B(ITR,1),LDBX)
ENDIF
ELSE
INZEE=1
DO I=2,IC-1
CALL FFTSTP(LOT,NFFT,N,LOT,N,Z(1,INZEE), &
Z(1,3-INZEE),TRIG,NOW(I),AFTER(I), &
BEFORE(I),ISIG)
INZEE=3-INZEE
ENDDO
IF(TRANSB == 'N'.OR.TRANSB == 'n') THEN
CALL FFTROT(LOT,NFFT,N,NFFT,LDBX,Z(1,INZEE), &
B(1,ITR),TRIG,NOW(IC),AFTER(IC),BEFORE(IC),ISIG)
ELSE
CALL FFTSTP(LOT,NFFT,N,LDBX,N,Z(1,INZEE), &
B(ITR,1),TRIG,NOW(IC),AFTER(IC),BEFORE(IC),ISIG)
ENDIF
ENDIF
ENDDO
IF(TRANSB == 'N'.OR.TRANSB == 'n') THEN
B(1:LDBX,M+1:LDBY) = CMPLX(0._dbl,0._dbl,dbl)
B(N+1:LDBX,1:M) = CMPLX(0._dbl,0._dbl,dbl)
ELSE
B(1:LDBX,N+1:LDBY) = CMPLX(0._dbl,0._dbl,dbl)
B(M+1:LDBX,1:M) = CMPLX(0._dbl,0._dbl,dbl)
ENDIF
!******************************************************************************
CONTAINS
!******************************************************************************
SUBROUTINE matmov ( n, m, a, lda, b, ldb )
IMPLICIT NONE
INTEGER :: n, m, lda, ldb
COMPLEX (dbl) :: a ( lda, * ), b ( ldb, * )
b ( 1:n , 1:m ) = a ( 1:n, 1:m )
END SUBROUTINE matmov
SUBROUTINE zgetmo ( a, lda, m, n, b, ldb )
IMPLICIT NONE
INTEGER :: lda, m, n, ldb
COMPLEX(dbl) :: a ( lda, n ), b ( ldb, m )
b ( 1:n, 1:m ) = TRANSPOSE ( a ( 1:m, 1:n ) )
END SUBROUTINE zgetmo
FUNCTION get_cache_size ( ) RESULT ( ncache )
IMPLICIT NONE
INTEGER ncache
#if defined ( __T3E )
ncache = 1024*8
#elif defined ( __SX5 ) || defined ( __T90 )
ncache = 1024*128
#elif defined ( __ALPHA )
ncache = 1024*8
#elif defined ( __SGI )
ncache = 1024*4
#elif defined ( __POWER2 )
ncache = 1024*10
#elif defined ( __HP )
ncache = 1024*64
#else
ncache = 1024*2
#endif
END FUNCTION get_cache_size
!******************************************************************************
END SUBROUTINE mltfftsg

352
src/library_tests.F Normal file
View file

@ -0,0 +1,352 @@
!-----------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
MODULE library_tests
USE kinds, ONLY : dbl
USE fft_tools, ONLY : init_fft, get_fft_library, fft3d, &
fft_radix_operations, FWFFT, BWFFT, FFT_RADIX_CLOSEST
USE global_types, ONLY : global_environment_type
USE machine, ONLY : m_cputime
USE stop_program, ONLY : stop_memory
IMPLICIT NONE
PRIVATE
PUBLIC :: lib_test
INTEGER :: runtest ( 100 )
!******************************************************************************
CONTAINS
!******************************************************************************
SUBROUTINE lib_test ( globenv )
TYPE ( global_environment_type ), INTENT ( IN ) :: globenv
INTEGER :: iw
iw = globenv % scr
IF ( globenv % ionode ) THEN
WRITE ( iw, '(T2,79("*"))' )
WRITE ( iw, '(A,T31,A,T80,A)' ) ' *',' PERFORMANCE TESTS ','*'
WRITE ( iw, '(T2,79("*"))' )
END IF
!
! CALL test_input ( globenv )
runtest = 0
runtest ( 3 ) = 1
!
IF ( runtest ( 1 ) == 1 ) CALL copy_test ( globenv )
!
IF ( runtest ( 2 ) == 1 ) CALL matmul_test ( globenv )
!
IF ( runtest ( 3 ) == 1 ) CALL fft_test ( globenv )
END SUBROUTINE lib_test
!******************************************************************************
SUBROUTINE copy_test ( globenv )
TYPE ( global_environment_type ), INTENT ( IN ) :: globenv
REAL ( dbl ), DIMENSION ( : ), ALLOCATABLE :: ca, cb
INTEGER :: len, ntim, iw, ierr, i, j
REAL ( dbl ) :: perf, tstart, tend, t
! test for copy --> Cache size
iw = globenv % scr
IF ( globenv % ionode ) WRITE ( iw, '(//,A,/)' ) " Test of copy ( F95 ) "
DO i = 6, 24
len = 2**i
ALLOCATE ( ca ( len ), STAT = ierr )
IF ( ierr /= 0 ) EXIT
ALLOCATE ( cb ( len ), STAT = ierr )
IF ( ierr /= 0 ) EXIT
CALL random_number ( ca )
ntim = NINT ( 1.e7_dbl / REAL ( len, dbl ) )
ntim = MAX ( ntim, 1 )
ntim = MIN ( ntim, 50000 )
tstart = m_cputime ( )
DO j = 1, ntim
cb ( : ) = ca ( : )
ca ( 1 ) = REAL ( j, dbl )
END DO
tend = m_cputime ( )
t = tend - tstart
perf = REAL ( ntim, dbl ) * REAL ( len, dbl ) * 1.e-6_dbl / t
IF ( globenv % ionode ) THEN
WRITE ( iw, '(A,i2,i10,A,T59,F14.4,A)' ) " Copy test: Size = 2^",i, &
len/1024," Kwords",perf," Mcopy/s"
END IF
DEALLOCATE ( ca , STAT = ierr )
IF ( ierr /= 0 ) CALL stop_memory ( "lib_test", "ca" )
DEALLOCATE ( cb , STAT = ierr )
IF ( ierr /= 0 ) CALL stop_memory ( "lib_test", "cb" )
END DO
END SUBROUTINE copy_test
!******************************************************************************
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
REAL ( dbl ) :: perf, tstart, tend, t
! test for matrix multpies
iw = globenv % scr
IF ( globenv % ionode ) WRITE ( iw, '(//,A,/)' ) " Test of matmul ( F95 ) "
DO i = 4, 10, 2
len = 2**i + 1
ALLOCATE ( ma ( len, len ), STAT = ierr )
IF ( ierr /= 0 ) EXIT
ALLOCATE ( mb ( len, len ), STAT = ierr )
IF ( ierr /= 0 ) EXIT
ALLOCATE ( mc ( len, len ), STAT = ierr )
IF ( ierr /= 0 ) EXIT
CALL random_number ( ma )
CALL random_number ( mb )
ntim = NINT ( 1.e8_dbl / ( 2._dbl * REAL ( len, dbl )**3 ) )
ntim = MAX ( ntim, 1 )
ntim = MIN ( ntim, 1000 )
tstart = m_cputime ( )
DO j = 1, ntim
mc = matmul ( ma, mb )
ma ( 1, 1 ) = REAL ( j, dbl )
END DO
tend = m_cputime ( )
t = tend - tstart
perf = REAL ( ntim, dbl ) * 2._dbl * REAL ( len, dbl )**3 * 1.e-6_dbl / t
IF ( globenv % ionode ) THEN
WRITE ( iw, '(A,i6,T59,F14.4,A)' ) &
" Matrix multiply test: c = a * b Size = ",len, perf," Mflop/s"
END IF
tstart = m_cputime ( )
DO j = 1, ntim
ma = matmul ( ma, mb )
ma ( 1, 1 ) = REAL ( j, dbl )
END DO
tend = m_cputime ( )
t = tend - tstart
perf = REAL ( ntim, dbl ) * 2._dbl * REAL ( len, dbl )**3 * 1.e-6_dbl / t
IF ( globenv % ionode ) THEN
WRITE ( iw, '(A,i6,T59,F14.4,A)' ) &
" Matrix multiply test: a = a * b Size = ",len, perf," Mflop/s"
END IF
tstart = m_cputime ( )
DO j = 1, ntim
mc = matmul ( ma, transpose ( mb ) )
ma ( 1, 1 ) = REAL ( j, dbl )
END DO
tend = m_cputime ( )
t = tend - tstart
perf = REAL ( ntim, dbl ) * 2._dbl * REAL ( len, dbl )**3 * 1.e-6_dbl / t
IF ( globenv % ionode ) THEN
WRITE ( iw, '(A,i6,T59,F14.4,A)' ) &
" Matrix multiply test: c = a * b(T) Size = ",len, perf," Mflop/s"
END IF
tstart = m_cputime ( )
DO j = 1, ntim
mc = matmul ( transpose ( ma ), mb )
ma ( 1, 1 ) = REAL ( j, dbl )
END DO
tend = m_cputime ( )
t = tend - tstart
perf = REAL ( ntim, dbl ) * 2._dbl * REAL ( len, dbl )**3 * 1.e-6_dbl / t
IF ( globenv % ionode ) THEN
WRITE ( iw, '(A,i6,T59,F14.4,A)' ) &
" Matrix multiply test: c = a(T) * b Size = ",len, perf," Mflop/s"
END IF
DEALLOCATE ( ma , STAT = ierr )
IF ( ierr /= 0 ) CALL stop_memory ( "lib_test", "ma" )
DEALLOCATE ( mb , STAT = ierr )
IF ( ierr /= 0 ) CALL stop_memory ( "lib_test", "mb" )
DEALLOCATE ( mc , STAT = ierr )
IF ( ierr /= 0 ) CALL stop_memory ( "lib_test", "mc" )
END DO
! test for matrix multpies
IF ( globenv % ionode ) WRITE ( iw, '(//,A,/)' ) " Test of matmul ( BLAS ) "
DO i = 4, 10, 2
len = 2**i + 1
ALLOCATE ( ma ( len, len ), STAT = ierr )
IF ( ierr /= 0 ) EXIT
ALLOCATE ( mb ( len, len ), STAT = ierr )
IF ( ierr /= 0 ) EXIT
ALLOCATE ( mc ( len, len ), STAT = ierr )
IF ( ierr /= 0 ) EXIT
CALL random_number ( ma )
CALL random_number ( mb )
ntim = NINT ( 1.e8_dbl / ( 2._dbl * REAL ( len, dbl )**3 ) )
ntim = MAX ( ntim, 1 )
ntim = MIN ( ntim, 1000 )
tstart = m_cputime ( )
DO j = 1, ntim
CALL DGEMM ( "N", "N", len, len, len, 1._dbl, ma, len, mb, len, 0._dbl, mc, len )
END DO
tend = m_cputime ( )
t = tend - tstart
perf = REAL ( ntim, dbl ) * 2._dbl * REAL ( len, dbl )**3 * 1.e-6_dbl / t
IF ( globenv % ionode ) THEN
WRITE ( iw, '(A,i6,T59,F14.4,A)' ) &
" Matrix multiply test: c = a * b Size = ",len, perf," Mflop/s"
END IF
tstart = m_cputime ( )
DO j = 1, ntim
CALL DGEMM ( "N", "N", len, len, len, 1._dbl, ma, len, mb, len, 0._dbl, mc, len )
CALL DCOPY ( len * len , mc, 1, ma, 1)
END DO
tend = m_cputime ( )
t = tend - tstart
perf = REAL ( ntim, dbl ) * 2._dbl * REAL ( len, dbl )**3 * 1.e-6_dbl / t
IF ( globenv % ionode ) THEN
WRITE ( iw, '(A,i6,T59,F14.4,A)' ) &
" Matrix multiply test: a = a * b Size = ",len, perf," Mflop/s"
END IF
tstart = m_cputime ( )
DO j = 1, ntim
CALL DGEMM ( "N", "T", len, len, len, 1._dbl, ma, len, mb, len, 0._dbl, mc, len )
END DO
tend = m_cputime ( )
t = tend - tstart
perf = REAL ( ntim, dbl ) * 2._dbl * REAL ( len, dbl )**3 * 1.e-6_dbl / t
IF ( globenv % ionode ) THEN
WRITE ( iw, '(A,i6,T59,F14.4,A)' ) &
" Matrix multiply test: c = a * b(T) Size = ",len, perf," Mflop/s"
END IF
tstart = m_cputime ( )
DO j = 1, ntim
CALL DGEMM ( "T", "N", len, len, len, 1._dbl, ma, len, mb, len, 0._dbl, mc, len )
END DO
tend = m_cputime ( )
t = tend - tstart
perf = REAL ( ntim, dbl ) * 2._dbl * REAL ( len, dbl )**3 * 1.e-6_dbl / t
IF ( globenv % ionode ) THEN
WRITE ( iw, '(A,i6,T59,F14.4,A)' ) &
" Matrix multiply test: c = a(T) * b Size = ",len, perf," Mflop/s"
END IF
DEALLOCATE ( ma , STAT = ierr )
IF ( ierr /= 0 ) CALL stop_memory ( "lib_test", "ma" )
DEALLOCATE ( mb , STAT = ierr )
IF ( ierr /= 0 ) CALL stop_memory ( "lib_test", "mb" )
DEALLOCATE ( mc , STAT = ierr )
IF ( ierr /= 0 ) CALL stop_memory ( "lib_test", "mc" )
END DO
END SUBROUTINE matmul_test
!******************************************************************************
SUBROUTINE fft_test ( globenv )
TYPE ( global_environment_type ), INTENT ( IN ) :: globenv
COMPLEX ( dbl ), DIMENSION ( :, :, : ), ALLOCATABLE :: ca, cb
REAL ( dbl ), DIMENSION ( :, :, : ), ALLOCATABLE :: ra
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
! test for 3d FFT
iw = globenv % scr
IF ( globenv % ionode ) WRITE ( iw, '(//,A,/)' ) " Test of 3D-FFT "
DO iall = 1, 100
SELECT CASE ( iall )
CASE DEFAULT
EXIT
CASE ( 1 )
CALL init_fft ( "FFTSG" )
method = "FFTSG "
CASE ( 2 )
CALL init_fft ( "FFTW" )
method = "FFTW "
END SELECT
n = 4
CALL fft3d ( 1, n, zz, status=stat )
IF ( stat == 0 ) THEN
DO it = 1, 3
radix_in = ndate ( it )
CALL fft_radix_operations ( radix_in, radix_out, FFT_RADIX_CLOSEST )
len = radix_out
n = len
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 )
IF ( ierr == 0 ) THEN
CALL random_number ( ra )
ca = ra
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 = MAX ( ntim, 1 )
ntim = MIN ( ntim, 1000 )
scale = 1._dbl / REAL ( len**3, dbl )
tstart = m_cputime ( )
DO j = 1, ntim
CALL fft3d ( FWFFT, n, ca )
CALL fft3d ( BWFFT, n, ca, SCALE = scale )
END DO
tend = m_cputime ( )
t = tend - tstart
perf = REAL ( ntim, dbl ) * 2._dbl * flops * 1.e-6_dbl / t
IF ( globenv % ionode ) THEN
WRITE ( iw, '(T2,A,A,i6,T59,F14.4,A)' ) &
adjustr(method)," test (in-place) Size = ",len, perf," Mflop/s"
END IF
DEALLOCATE ( ca , STAT = ierr )
IF ( ierr /= 0 ) CALL stop_memory ( "lib_test", "ca" )
DEALLOCATE ( ra , STAT = ierr )
IF ( ierr /= 0 ) CALL stop_memory ( "lib_test", "ra" )
END IF
END DO
END IF
IF ( globenv % ionode ) WRITE ( iw, * )
END DO
END SUBROUTINE fft_test
!******************************************************************************
END MODULE library_tests
!******************************************************************************

7
src/linklist_control.F
View file

@ -19,6 +19,7 @@ MODULE linklist_control
USE particle_types, ONLY : particle_type
USE simulation_cell, ONLY : cell_type
USE stop_program, ONLY : stop_prg, stop_memory
USE timings, ONLY : timeset, timestop
IMPLICIT NONE
@ -162,11 +163,13 @@ SUBROUTINE list_control ( pnode, part, box )
REAL ( dbl ), ALLOCATABLE, SAVE :: r_last_update ( :, : )
REAL ( dbl ) :: displace, max_displace
LOGICAL :: list_update_flag, ionode
INTEGER :: i, nnodes, isos, iw
INTEGER :: i, nnodes, isos, iw, handle
CHARACTER ( LEN = 36 ) :: string
!------------------------------------------------------------------------------
CALL timeset ( 'LIST_CONTROL','I',' ',handle )
nnodes = SIZE ( pnode )
IF ( .NOT. ALLOCATED ( r_last_update ) ) THEN
ALLOCATE ( r_last_update ( 3, nnodes ), STAT = isos )
@ -234,6 +237,8 @@ SUBROUTINE list_control ( pnode, part, box )
END IF
counter = counter + 1
CALL timestop ( 0._dbl, handle )
END SUBROUTINE list_control

1
src/linklists.F
View file

@ -135,6 +135,7 @@ SUBROUTINE exclusion(molecule,pnode)
bond(2) = ll_bond%p2%iatom
! checking to see if jatom is bonded
match=.FALSE.
DO ii = 1, 2
IF (jatom==bond(ii)) THEN
match = .TRUE.

12
src/machine.F
View file

@ -31,4 +31,16 @@ MODULE machine
PUBLIC :: m_walltime, m_cputime, m_datum
PUBLIC :: m_hostnm, m_getcwd, m_getlog, m_getuid, m_getpid, m_getarg
!
! FUNCTION m_cputime()
! FUNCTION m_walltime()
! SUBROUTINE m_datum(cal_date)
! SUBROUTINE m_hostnm(hname)
! SUBROUTINE m_getcwd(curdir)
! SUBROUTINE m_getlog(user)
! SUBROUTINE m_getuid(uid)
! SUBROUTINE m_getpid(pid)
! SUBROUTINE m_getarg(i,arg)
!
END MODULE machine

31
src/nose.F
View file

@ -799,6 +799,8 @@ SUBROUTINE yoshida_coef ( nhcp, dt )
!------------------------------------------------------------------------------
SELECT CASE (nhcp % nyosh)
CASE DEFAULT
CALL stop_prg ( 'yoshida_coef', 'Value not available' )
CASE (1)
yosh_wt(1) = 1.0_dbl
CASE (3)
@ -819,8 +821,35 @@ SUBROUTINE yoshida_coef ( nhcp, dt )
yosh_wt(5) = yosh_wt(3)
yosh_wt(6) = yosh_wt(2)
yosh_wt(7) = yosh_wt(1)
CASE (9)
yosh_wt(1) = 0.192_dbl
yosh_wt(2) = 0.554910818409783619692725006662999_dbl
yosh_wt(3) = 0.124659619941888644216504240951585_dbl
yosh_wt(4) = -0.843182063596933505315033808282941_dbl
yosh_wt(5) = 1.0_dbl - 2.0_dbl*(yosh_wt(1)+yosh_wt(2)+&
yosh_wt(3)+yosh_wt(4))
yosh_wt(6) = yosh_wt(4)
yosh_wt(7) = yosh_wt(3)
yosh_wt(8) = yosh_wt(2)
yosh_wt(9) = yosh_wt(1)
CASE (15)
yosh_wt(1) = 0.102799849391985_dbl
yosh_wt(2) = -0.196061023297549e1_dbl
yosh_wt(3) = 0.193813913762276e1_dbl
yosh_wt(4) = -0.158240635368243_dbl
yosh_wt(5) = -0.144485223686048e1_dbl
yosh_wt(6) = 0.253693336566229_dbl
yosh_wt(7) = 0.914844246229740_dbl
yosh_wt(8) = 1.0_dbl - 2.0_dbl*(yosh_wt(1)+yosh_wt(2)+&
yosh_wt(3)+yosh_wt(4)+yosh_wt(5)+yosh_wt(6)+yosh_wt(7))
yosh_wt(9) = yosh_wt(7)
yosh_wt(10) = yosh_wt(6)
yosh_wt(11) = yosh_wt(5)
yosh_wt(12) = yosh_wt(4)
yosh_wt(13) = yosh_wt(3)
yosh_wt(14) = yosh_wt(2)
yosh_wt(15) = yosh_wt(1)
END SELECT
nhcp % dt_yosh = dt * yosh_wt / REAL ( nhcp % nc, dbl )
END SUBROUTINE yoshida_coef

10
src/pw_grids.F
View file

@ -54,11 +54,11 @@ SUBROUTINE pw_grid_setup ( cell, pw_grid, cutoff )
pw_grid % ngpts = PRODUCT ( pw_grid % npts )
IF ( PRESENT ( cutoff ) ) THEN
write ( 6, '( A, I8 )' ) "pw_grid_setup: #gpts = ", pw_grid % ngpts_cut
ELSE
write ( 6, '( A, I8 )' ) "pw_grid_setup: #gpts = ", pw_grid % ngpts
END IF
! IF ( PRESENT ( cutoff ) ) THEN
! write ( 6, '( A, I8 )' ) " PW_grid_setup| #gpts = ", pw_grid % ngpts_cut
! ELSE
! write ( 6, '( A, I8 )' ) " PW_grid_setup| #gpts = ", pw_grid % ngpts
! END IF
! Allocate the grid fields
CALL pw_grid_allocate ( pw_grid )

2
src/pws.F
View file

@ -7,7 +7,7 @@ MODULE pws
USE coefficient_types, ONLY : coeff_type, coeff_allocate, coeff_deallocate, &
SQUARE, SQUAREROOT
USE fft_tools, ONLY : fft3d => fft_wrap, &
USE fft_tools, ONLY : fft3d, &
FWFFTu => FWFFT, BWFFTu => BWFFT, FWFFTw => FWFFT, BWFFTw => BWFFT
USE kinds, ONLY : dbl
USE mathconstants, ONLY : fourpi

20
src/structure_types.F
View file

@ -9,20 +9,30 @@ MODULE structure_types
USE molecule_types, ONLY : molecule_structure_type, particle_node_type
USE particle_types, ONLY : particle_type
USE simulation_cell, ONLY : cell_type
USE pair_potential, ONLY : potentialparm_type
USE tbmd_types, ONLY : tbatom_param_type, tb_hopping_type
IMPLICIT NONE
PRIVATE
PUBLIC :: structure_type
PUBLIC :: structure_type, interaction_type
PUBLIC :: cell_type, particle_type, particle_node_type
PUBLIC :: molecule_structure_type
PUBLIC :: potentialparm_type, tbatom_param_type, tb_hopping_type
TYPE structure_type
CHARACTER ( LEN = 30 ) :: name
TYPE (cell_type ) :: box
TYPE (particle_type ), POINTER :: part ( : )
TYPE (particle_node_type ), POINTER :: pnode ( : )
TYPE (molecule_structure_type ), POINTER :: molecule ( : )
TYPE ( cell_type ) :: box
TYPE ( particle_type ), POINTER :: part ( : )
TYPE ( particle_node_type ), POINTER :: pnode ( : )
TYPE ( molecule_structure_type ), POINTER :: molecule ( : )
END TYPE structure_type
TYPE interaction_type
TYPE ( potentialparm_type ), DIMENSION ( :, : ), POINTER :: potparm
TYPE ( tbatom_param_type ), DIMENSION ( : ), POINTER :: tbatom
TYPE ( tb_hopping_type ), DIMENSION ( :, : ), POINTER :: tbhop
END TYPE interaction_type
END MODULE structure_types

221
src/tbmd_debug.F
View file

@ -6,21 +6,21 @@
MODULE tbmd_debug
USE ewald_parameters_types, ONLY : ewald_parameters_type
USE fist_force_numer, ONLY : force_nonbond_numer, force_recip_numer, &
ptens_numer, pvg_numer, potential_g_numer, de_g_numer, &
energy_recip_numer
USE tbmd_force, ONLY : force_control, debug_variables_type
USE global_types, ONLY : global_environment_type
USE kinds, ONLY : dbl
USE linklists, ONLY : dist_constraints, g3x3_constraints
USE md, ONLY : thermodynamic_type
USE molecule_types, ONLY : molecule_structure_type, particle_node_type
USE particle_types, ONLY : particle_type
USE pair_potential, ONLY : potentialparm_type, ener_coul
USE tbmd_force, ONLY : force_control, debug_variables_type
USE fist_force_numer, ONLY : force_nonbond_numer, force_recip_numer, &
ptens_numer, pvg_numer, potential_g_numer, de_g_numer, &
energy_recip_numer
USE linklists, ONLY : dist_constraints, g3x3_constraints
USE pw_grid_types, ONLY : pw_grid_type
USE pw_grids, ONLY : pw_grid_setup
USE pw_grid_types, ONLY : pw_grid_type, HALFSPACE
USE pw_grids, ONLY : pw_find_cutoff, pw_grid_setup
USE simulation_cell, ONLY : cell_type
USE stop_program, ONLY : stop_memory
USE stop_program, ONLY : stop_memory, stop_prg
IMPLICIT NONE
@ -51,16 +51,16 @@ SUBROUTINE debug_control ( globenv, ewald_param, part, pnode, molecule, box, &
! Locals
TYPE ( debug_variables_type ) :: dbg
TYPE ( pw_grid_type ) :: pw_grid
REAL ( dbl ), DIMENSION ( :, : ), ALLOCATABLE :: f_numer
REAL ( dbl ), DIMENSION ( :, : ), ALLOCATABLE :: rel, diff
INTEGER :: iflag, i, natoms, isos, iatom, iw, ir
TYPE ( pw_grid_type ) :: ewald_grid
INTEGER :: iflag, i, natoms, isos, iatom, iw, ir, npts_s(3), gmax
INTEGER, DIMENSION ( 2 ) :: dum
REAL ( dbl ) :: delta, energy_numer
REAL ( dbl ) :: delta, energy_numer, cutoff
REAL ( dbl ) :: e_numer, pv_test ( 3, 3 )
REAL ( dbl ) :: err1, numer, denom1, vec ( 3 ), e_bc, e_real, energy_tot
REAL ( dbl ) :: denom2, err2
REAL ( dbl ), DIMENSION ( :, : ), ALLOCATABLE :: f_bc, f_real
REAL ( dbl ), DIMENSION ( :, : ), ALLOCATABLE :: f_numer
REAL ( dbl ), DIMENSION ( :, : ), ALLOCATABLE :: rel, diff
!------------------------------------------------------------------------------
@ -107,8 +107,8 @@ SUBROUTINE debug_control ( globenv, ewald_param, part, pnode, molecule, box, &
ELSE
WRITE ( iw, '( A, T71, F10.6 )' ) ' DELTA = ', delta
END IF
CALL force_nonbond_numer ( pnode, box, potparm, delta, f_numer, &
energy_numer )
CALL force_nonbond_numer ( ewald_param, pnode, box, potparm, &
delta, f_numer, energy_numer )
WRITE ( iw, '( A, T61, E20.14 )' ) ' NON BOND NUMER ENERGY = ', &
energy_numer
WRITE ( iw, '( A, T61, E20.14 )' ) ' NON BOND ANAL ENERGY = ', &
@ -144,101 +144,112 @@ SUBROUTINE debug_control ( globenv, ewald_param, part, pnode, molecule, box, &
END IF
! Debug g-space
! IF ( ewald_param % ewald_type /= 'NONE' ) THEN
! WRITE ( iw, '( A )' ) ' DO YOU WANT TO DEBUG YOUR G-SPACE (1=yes)?'
! READ ( ir, * ) iflag
! IF ( iflag == 1 ) THEN
! CALL pw_grid_setup ( box, pw_grid, ewald_param % gmax )
! WRITE ( iw, '( A )' ) &
! ' 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= ', &
! pw_grid % ngpts_cut
! WRITE ( iw, '( A, T71, F10.4 )' ) 'ALPHA= ', &
! ewald_param % alpha
! CALL energy_recip_numer ( pnode, box, potparm, pw_grid, &
! energy_numer, pw_grid % ngpts_cut )
! WRITE ( iw, '( A, T61, G20.14 )' ) 'G-SPACE ANAL ENERGY = ', &
! dbg % pot_g
! WRITE ( iw, '( A, T61, G20.14 )' ) &
! 'G-SPACE NUMERICAL ENERGY = ', energy_numer
! END IF
! WRITE ( iw, '( A )' ) &
! ' DO YOU WANT TO DEBUG YOUR G-SPACE FORCES (1=yes)?'
! READ ( ir, * ) iflag
! IF ( iflag == 1 ) THEN
! 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
! WRITE ( iw, '( A, T71, F10.6 )' ) &
! ' DELTA (changed to default) = ', delta
! ELSE
! WRITE ( iw, '( A, T71, F10.6 )' ) ' DELTA = ', delta
! END IF
!
! CALL force_recip_numer ( pnode, box, potparm, gvec, delta, f_numer )
!
IF ( ewald_param % ewald_type /= 'NONE' ) THEN
WRITE ( iw, '( A )' ) ' DO YOU WANT TO DEBUG YOUR G-SPACE (1=yes)?'
READ ( ir, * ) iflag
IF ( iflag == 1 ) THEN
gmax = ewald_param % gmax
ewald_grid % bounds ( 1, : ) = -gmax / 2
ewald_grid % bounds ( 2, : ) = +gmax / 2
npts_s = (/ gmax, gmax, gmax /)
ewald_grid % grid_span = HALFSPACE
CALL pw_find_cutoff ( npts_s, box, cutoff )
CALL pw_grid_setup( box, ewald_grid, cutoff)
WRITE ( iw, '( A )' ) &
' 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= ', &
ewald_grid % ngpts_cut
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 = ', &
dbg % pot_g
WRITE ( iw, '( A, T61, G20.14 )' ) &
'G-SPACE NUMERICAL ENERGY = ', energy_numer
END IF
WRITE ( iw, '( A )' ) &
' DO YOU WANT TO DEBUG YOUR G-SPACE FORCES (1=yes)?'
READ ( ir, * ) iflag
IF ( iflag == 1 ) THEN
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
WRITE ( iw, '( A, T71, F10.6 )' ) &
' DELTA (changed to default) = ', delta
ELSE
WRITE ( iw, '( A, T71, F10.6 )' ) ' DELTA = ', delta
END IF
CALL force_recip_numer ( ewald_param, pnode, box, ewald_grid, &
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
!
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 )' ) &
! 'MAXIMUM ABSOLUTE ERROR = ', maxval(diff)
! WRITE ( iw, '( A,T61,G20.14 )' ) &
! 'MAXIMUM RELATIVE ERROR = ', maxval(rel)
WRITE ( iw, '( A,T61,G20.14 )' ) &
'MAXIMUM ABSOLUTE ERROR = ', maxval(diff)
WRITE ( iw, '( A,T61,G20.14 )' ) &
'MAXIMUM RELATIVE ERROR = ', maxval(rel)
!
! write out the particle number and forces of
! the max absolute and relative error
! 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 =', &
! dbg % f_g(:,dum(2))
! WRITE ( iw, '( A,T21,3G20.14 )' ) ' 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 =', &
! dbg % f_g(:,dum(2))
! WRITE ( iw, '( A,T21,3G20.14 )' ) ' F_NUMR G-SPACE =', &
! f_numer(:,dum(2))
! END IF
!
! WRITE ( iw, '( A )' ) &
! ' DO YOU WANT TO DEBUG YOUR G-SPACE VIRIAL (1=yes)?'
! READ ( ir, * ) iflag
! IF ( iflag == 1 ) THEN
! write out the particle number and forces of the max absolute
! and relative error
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 =', &
dbg % f_g(:,dum(2))
WRITE ( iw, '( A,T21,3G20.14 )' ) ' 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 =', &
dbg % f_g(:,dum(2))
WRITE ( iw, '( A,T21,3G20.14 )' ) ' F_NUMR G-SPACE =', &
f_numer(:,dum(2))
END IF
WRITE ( iw, '( A )' ) &
' DO YOU WANT TO DEBUG YOUR G-SPACE VIRIAL (1=yes)?'
READ ( ir, * ) iflag
IF ( iflag == 1 ) THEN
!
! get numerical virial
! 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
! WRITE ( iw, '( A,T71,F10.6 )' ) &
! ' DELTA (changed to default) = ', delta
! ELSE
! WRITE ( iw, '( A,T71,F10.6 )' ) ' DELTA = ', delta
! END IF
! CALL pvg_numer ( pnode, box, potparm, gvec, pv_test, delta )
!
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
WRITE ( iw, '( A,T71,F10.6 )' ) &
' DELTA (changed to default) = ', delta
ELSE
WRITE ( iw, '( A,T71,F10.6 )' ) ' DELTA = ', delta
END IF
CALL pvg_numer ( ewald_param, pnode, box, &
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,:)
! END DO
! WRITE ( iw, '( A )' ) ' PV G-SPACE'
! DO i = 1, 3
! WRITE ( iw, '( T21,3G20.14 )' ) dbg % pv_g(i,:)
! END DO
! END IF
! END IF
! END IF
!
WRITE ( iw, '( A )' ) ' PV G-SPACE NUMERICAL'
DO i = 1, 3
WRITE ( iw, '( T21,3G20.14 )' ) pv_test(i,:)
END DO
WRITE ( iw, '( A )' ) ' PV G-SPACE'
DO i = 1, 3
WRITE ( iw, '( T21,3G20.14 )' ) dbg % pv_g(i,:)
END DO
END IF
END IF
END IF
! Debug real-space virial
WRITE ( iw, '( A )' ) &
' DO YOU WANT TO DEBUG YOUR NON BOND (REAL SPACE) VIRIAL (1=yes)?'
@ -255,7 +266,7 @@ SUBROUTINE debug_control ( globenv, ewald_param, part, pnode, molecule, box, &
ELSE
WRITE ( iw, '( A, T71, F10.6 )' ) ' DELTA = ', delta
END IF
CALL ptens_numer ( pnode, box, potparm, pv_test, delta )
CALL ptens_numer ( ewald_param, pnode, box, potparm, pv_test, delta )
!
! writing out the virials
WRITE ( iw, '( A )' ) ' PV NUMERICAL'

9
src/tbmd_force.F
View file

@ -87,12 +87,11 @@ SUBROUTINE force_control ( molecule, pnode, part, box, thermo, &
IF ( first_time ) THEN
SELECT CASE ( ewald_param % ewald_type )
CASE ( 'EWALD_GAUSS' )
CALL ewald_initialize ( dg, part, pnode, fc_global % group, &
ewald_param, box, thermo, fc_global % scr, &
ewald_grid = grid_ewald )
CALL ewald_initialize ( dg, part, pnode, fc_global, &
ewald_param, box, thermo, ewald_grid = grid_ewald )
CASE ( 'PME_GAUSS' )
CALL ewald_initialize ( dg, part, pnode, fc_global % group, &
ewald_param, box, thermo, fc_global % scr, &
CALL ewald_initialize ( dg, part, pnode, fc_global, &
ewald_param, box, thermo, &
pme_small_grid = grid_s, pme_big_grid = grid_b )
END SELECT
END IF

129
src/tbmd_initialize.F Normal file
View file

@ -0,0 +1,129 @@
!------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!------------------------------------------------------------------------------!
!
MODULE tbmd_initialize
!
!------------------------------------------------------------------------------!
USE kinds, ONLY : dbl
USE stop_program, ONLY : stop_memory
USE tbmd_types, ONLY : tbatom_param_type, tb_hopping_type
USE particle_types, ONLY : particle_type
!
IMPLICIT NONE
!
PRIVATE
PUBLIC :: tbmd_init, tb_get_numel, tb_get_basize
!
!------------------------------------------------------------------------------!
!
CONTAINS
!
!------------------------------------------------------------------------------!
SUBROUTINE tbmd_init(part,tbatom,tbhop)
IMPLICIT NONE
TYPE (particle_type), INTENT (IN) :: part(:)
TYPE (tbatom_param_type), DIMENSION(:), INTENT (INOUT) :: tbatom
TYPE (tb_hopping_type), DIMENSION(:,:), INTENT (INOUT) :: tbhop
INTEGER :: iat, jat, i, j, k, n, ii, ni, nj, is, js, ju, ios
INTEGER :: natom_type
!
! count number of atoms per atom type
natom_type = size(tbatom)
DO iat = 1, natom_type
tbatom(iat) %natoms = 0
END DO
DO iat = 1, size(part)
i = part(iat) %prop%ptype
tbatom(i) %natoms = tbatom(i) %natoms + 1
END DO
! count number of basis function per atom type
! set up diagonal part of hamiltonian
DO iat = 1, natom_type
n = tbatom(iat) %bsf(0) + 3*tbatom(iat) %bsf(1) + &
5*tbatom(iat) %bsf(2) + 7*tbatom(iat) %bsf(3)
tbatom(iat) %nbsf = n
ALLOCATE (tbatom(iat)%h_diag(n),STAT=ios)
IF (ios/=0) CALL stop_memory('tbmd_initialize', &
'tbatom(iat)%h_diag',n)
ii = 0
DO i = 0, 3
DO j = 1, tbatom(iat) %bsf(i)
tbatom(iat) %h_diag(ii+1:ii+2*i+1) = tbatom(iat) &
%orbital_energy(i,j)
ii = ii + 2*i + 1
END DO
END DO
END DO
! set up shell pointers
DO iat = 1, natom_type
n = tbatom(iat) %bsf(0) + tbatom(iat) %bsf(1) + &
tbatom(iat) %bsf(2) + tbatom(iat) %bsf(3)
ALLOCATE (tbatom(iat)%spoint(n,2),STAT=ios)
IF (ios/=0) CALL stop_memory('tbmd_initialize', &
'tbatom(iat)%spoint',2*n)
ii = 0
DO i = 0, 3
DO j = 1, tbatom(iat) %bsf(i)
ii = ii + 1
tbatom(iat) %spoint(ii,2) = i
END DO
END DO
tbatom(iat) %spoint(1,1) = 0
DO j = 2, n
ii = tbatom(iat) %spoint(j-1,2)
tbatom(iat) %spoint(j,1) = tbatom(iat) %spoint(j-1,1) + 2*ii + 1
END DO
END DO
! pointer: shell pair to hopping potential
DO iat = 1, natom_type
DO jat = 1, natom_type
ni = size(tbatom(iat)%spoint(:,1))
nj = size(tbatom(jat)%spoint(:,1))
ALLOCATE (tbhop(iat,jat)%llpoint(ni,nj),STAT=ios)
IF (ios/=0) CALL stop_memory('tbmd_initialize', &
'tbatom(iat)%llpoint',ni*nj)
k = 0
DO i = 1, ni
is = tbatom(iat) %spoint(i,2)
ju = 1
IF (iat==jat) ju = i
DO j = ju, nj
js = tbatom(jat) %spoint(j,2)
tbhop(iat,jat) %llpoint(i,j) = k + 1
k = k + min(is,js)
END DO
END DO
END DO
END DO
!
END SUBROUTINE tbmd_init
!------------------------------------------------------------------------------!
FUNCTION tb_get_basize(tbatom) RESULT (basis_size)
IMPLICIT NONE
TYPE (tbatom_param_type), DIMENSION(:), INTENT (IN) :: tbatom
INTEGER :: basis_size, i
basis_size = 0
DO i = 1, size(tbatom)
basis_size = basis_size + tbatom(i) %nbsf*tbatom(i) %natoms
END DO
END FUNCTION tb_get_basize
!------------------------------------------------------------------------------!
FUNCTION tb_get_numel(tbatom,charge) RESULT (number_of_electrons)
IMPLICIT NONE
TYPE (tbatom_param_type), DIMENSION(:), INTENT (IN) :: tbatom
REAL (dbl) , INTENT(IN) :: charge
INTEGER :: number_of_electrons, i
number_of_electrons = 0
DO i = 1, size(tbatom)
number_of_electrons = number_of_electrons + &
tbatom(i) %valence_charge*tbatom(i) %natoms
END DO
number_of_electrons = number_of_electrons + charge
END FUNCTION tb_get_numel
!------------------------------------------------------------------------------!
END MODULE tbmd_initialize
!------------------------------------------------------------------------------!

123
src/tbmd_input.F
View file

@ -5,138 +5,27 @@
MODULE tbmd_input
!------------------------------------------------------------------------------!
USE kinds, ONLY : dbl
USE tbmd_global, ONLY : tbatom, tbhop
USE tbmd_types, ONLY : tbatom_param_type, tb_hopping_type
USE stop_program, ONLY : stop_prg, stop_memory
USE util, ONLY : get_unit
USE string_utilities, ONLY : uppercase, xstring, &
str_search, str_comp, make_tuple
USE parser, ONLY : parser_init, parser_end, read_line, test_next, &
cfield, p_error, get_real, get_int, stop_parser
USE message_passing, ONLY : mp_bcast
USE ewald_parameters_types, ONLY : ewald_parameters_type
USE input_types, ONLY : setup_parameters_type
USE global_types, ONLY : global_environment_type
IMPLICIT NONE
!
PRIVATE
PUBLIC :: read_tbmd_section, read_tb_hamiltonian, &
read_tb_hopping_elements
PUBLIC :: read_tb_hamiltonian, read_tb_hopping_elements
!
!------------------------------------------------------------------------------!
!
CONTAINS
!
!!>----------------------------------------------------------------------------!
!! SECTION: &tbmd ... &end !
!! !
!! simulation [energy,md,debug] !
!! units [atomic] !
!! periodic [0,1][0,1][0,1] !
!! set_file "filename" !
!! input_file "filename" !
!! !
!!tmptmptmptmptmp (to be included)????
!! symmetry [on,off] !
!!<----------------------------------------------------------------------------!
SUBROUTINE read_tbmd_section ( setup, ewald_param, tbmdpar )
IMPLICIT NONE
! Arguments
TYPE(setup_parameters_type), INTENT(INOUT) :: setup
TYPE(ewald_parameters_type), INTENT(INOUT) :: ewald_param
TYPE(global_environment_type), INTENT(INOUT) :: tbmdpar
! Locals
INTEGER :: ierror, ilen, ia, ie, i, j, n, iw, source, group
CHARACTER (len=20) :: string
CHARACTER (len=5) :: label, str2
CHARACTER (len=3), PARAMETER :: yn(0:1) = (/ ' NO', 'YES'/)
!------------------------------------------------------------------------------!
!..defaults
setup % run_type = 'ENERGY'
setup % unit_type = 'ATOMIC'
setup % perd = 1
ewald_param % alpha = 0.4_dbl
ewald_param % gmax = 10
ewald_param % ns_max = 10
ewald_param % epsilon = 1.e-6_dbl
ewald_param % ewald_type = 'NONE'
CALL xstring(tbmdpar%project_name,ia,ie)
setup % set_file_name = tbmdpar%project_name(ia:ie) // '.set'
setup % input_file_name = tbmdpar%project_name(ia:ie) // '.dat'
iw = tbmdpar % scr
!..parse the input section
label = '&TBMD'
CALL parser_init(tbmdpar%input_file_name,label,ierror,tbmdpar)
IF (ierror/=0) THEN
IF (tbmdpar%ionode) &
WRITE (iw,'(a)') ' No input section &TBMD found '
ELSE
CALL read_line
DO WHILE (test_next()/='X')
ilen = 8
CALL cfield(string,ilen)
CALL uppercase(string)
SELECT CASE (string)
CASE DEFAULT
CALL p_error()
CALL stop_parser('read_tbmd','unknown option')
CASE ('SIMULATI')
ilen = 20
CALL cfield(setup%run_type,ilen)
CALL uppercase(setup%run_type)
CASE ('PRINTLEV')
tbmdpar%print_level = get_int()
CASE ('UNITS')
ilen = 20
CALL cfield(setup%unit_type,ilen)
CALL uppercase(setup%unit_type)
CASE ('PERIODIC')
setup%perd(1) = get_int()
setup%perd(2) = get_int()
setup%perd(3) = get_int()
CASE ('SET_FILE')
ilen = 20
CALL cfield(setup%set_file_name,ilen)
CASE ('INP_FILE')
ilen = 20
CALL cfield(setup%input_file_name,ilen)
END SELECT
CALL read_line
END DO
END IF
CALL parser_end
!..end of parsing the input section
!..write some information to output
IF (tbmdpar%ionode) THEN
IF (tbmdpar%print_level>=0) THEN
WRITE (iw,'(A,T71,A)') ' TBMD| Run type ', adjustr(setup%run_type)
WRITE (iw,'(A,T71,A)') ' TBMD| Unit type ', adjustr(setup%unit_type)
WRITE (iw,'(A,T78,A)') ' TBMD| Periodic in X direction ', &
yn(setup%perd(1))
WRITE (iw,'(A,T78,A)') ' TBMD| Periodic in Y direction ', &
yn(setup%perd(2))
WRITE (iw,'(A,T78,A)') ' TBMD| Periodic in Z direction ', &
yn(setup%perd(3))
WRITE (iw,'(A,T61,A)') ' TBMD| Set file name', &
adjustr(setup%set_file_name)
WRITE (iw,'(A,T61,A)') ' TBMD| Input file name', &
adjustr(setup%input_file_name)
WRITE (iw,'(A,T76,I5)') ' TBMD| Print level ', tbmdpar%print_level
WRITE (iw,*)
END IF
END IF
END SUBROUTINE read_tbmd_section
!------------------------------------------------------------------------------!
SUBROUTINE read_tb_hamiltonian ( setup, atom_names, tbmdpar )
SUBROUTINE read_tb_hamiltonian ( setup, atom_names, tbatom, tbmdpar )
IMPLICIT NONE
@ -144,6 +33,8 @@ SUBROUTINE read_tb_hamiltonian ( setup, atom_names, tbmdpar )
TYPE(setup_parameters_type), INTENT(INOUT) :: setup
TYPE(global_environment_type), INTENT(IN) :: tbmdpar
CHARACTER ( LEN = * ), INTENT(IN) :: atom_names ( : )
TYPE (tbatom_param_type), DIMENSION (:), POINTER :: tbatom
! Locals
INTEGER :: ierror, ilen, i, j, k, n, iw, source, group, iat, nt, ios, i1
@ -263,7 +154,7 @@ SUBROUTINE set_basis(basis,bsf)
END SUBROUTINE set_basis
!------------------------------------------------------------------------------!
SUBROUTINE read_tb_hopping_elements ( setup, atom_names, tbmdpar )
SUBROUTINE read_tb_hopping_elements ( setup, atom_names, tbatom, tbhop, tbmdpar )
IMPLICIT NONE
@ -271,6 +162,8 @@ SUBROUTINE read_tb_hopping_elements ( setup, atom_names, tbmdpar )
TYPE(setup_parameters_type), INTENT(INOUT) :: setup
TYPE(global_environment_type), INTENT(INOUT) :: tbmdpar
CHARACTER ( LEN = * ), INTENT(IN) :: atom_names ( : )
TYPE (tbatom_param_type), DIMENSION (:) :: tbatom
TYPE (tb_hopping_type), DIMENSION (:,:), POINTER :: tbhop
INTEGER :: ierror, ilen, i, j, k, l, n, iw, source, group, iat, ios
INTEGER :: l1, l2, k1, k2, kk, nk, lu, ku