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

*** empty log message ***

Browse source 
svn-origin-rev: 33
This commit is contained in:
Gloria Tabacchi 2001-07-01 17:44:50 +00:00
parent 8b3557af19
commit a3094f855c
4 changed files with 1185 additions and 0 deletions
Download patch file Download diff file Expand all files Collapse all files

3
src/OBJECTDEFS
View file

@ -1,5 +1,6 @@
OBJECTS_GENERIC =\
ai_coulomb.o\
ai_derivatives.o\
ai_kinetic.o\
ai_nuclear.o\
ai_overlap.o\
@ -39,6 +40,7 @@ OBJECTS_GENERIC =\
diis.o\
dump.o\
eigenvalueproblems.o\
empirical_parameters.o\
environment.o\
erf_fn.o\
ewald_parameters_types.o\
@ -133,6 +135,7 @@ OBJECTS_GENERIC =\
pol_electrostatics.o\
pol_force.o\
pol_grids.o\
pol_overlap.o\
pol_setup.o\
potential_types.o\
potentials.o\

176
src/ai_derivatives.F Normal file
View file

@ -0,0 +1,176 @@
MODULE ai_derivatives
! Purpose: Calculate the first derivative of an integral block.
! Literature: S. Obara and A. Saika, J. Chem. Phys. 84, 3963 (1986)
! History: - Creation (05.10.2000, Matthias Krack)
! *****************************************************************************
! ax,ay,az : Angular momentum index numbers of orbital a.
! bx,by,bz : Angular momentum index numbers of orbital b.
! coset : Cartesian orbital set pointer.
! l{a,b} : Angular momentum quantum number of shell a or b.
! l{a,b}_max: Maximum angular momentum quantum number of shell a or b.
! l{a,b}_min: Minimum angular momentum quantum number of shell a or b.
! ncoset : Number of orbitals in a Cartesian orbital set.
! npgf{a,b} : Degree of contraction of shell a or b.
! rab : Distance vector between the atomic centers a and b.
! rab2 : Square of the distance between the atomic centers a and b.
! rac : Distance vector between the atomic centers a and c.
! rac2 : Square of the distance between the atomic centers a and c.
! rbc : Distance vector between the atomic centers b and c.
! rbc2 : Square of the distance between the atomic centers b and c.
! rpgf{a,b} : Radius of the primitive Gaussian-type function a or b.
! zet{a,b} : Exponents of the Gaussian-type functions a or b.
! zetp : Reciprocal of the sum of the exponents of orbital a and b.
! *****************************************************************************
USE kinds, ONLY: wp => dp
IMPLICIT NONE
PRIVATE
! *** Public subroutines ***
PUBLIC :: dabdr
! *****************************************************************************
CONTAINS
! *****************************************************************************
SUBROUTINE dabdr(la_max,zeta,rpgfa,la_min,&
lb_max,zetb,rpgfb,lb_min,&
dab,ab,dabdx,dabdy,dabdz)
! Purpose: Calculate the first derivative of an integral block.
! History: - Creation (05.10.2000, Matthias Krack)
! ***************************************************************************
USE orbital_pointers, ONLY: coset,ncoset
REAL(wp), INTENT(IN) :: dab
INTEGER, INTENT(IN) :: la_max,la_min,lb_max,lb_min
REAL(wp), DIMENSION(:), INTENT(IN) :: rpgfa,rpgfb,zeta,zetb
REAL(wp), DIMENSION(:,:), INTENT(IN) :: ab
REAL(wp), DIMENSION(:,:), INTENT(OUT) :: dabdx,dabdy,dabdz
! *** Local variables ***
REAL(wp) :: fa,fx,fy,fz
INTEGER :: ax,ay,az,bx,by,bz,coa,coamx,coamy,coamz,coapx,coapy,coapz,cob,&
coda,i,ipgf,j,jpgf,la,lb,na,nb,nda
! ---------------------------------------------------------------------------
! *** Loop over all pairs of primitive Gaussian-type functions ***
na = 0
nda = 0
DO ipgf=1,SIZE(zeta)
fa = 2.0_wp*zeta(ipgf)
nb = 0
DO jpgf=1,SIZE(zetb)
! *** Screening ***
IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) THEN
DO j=nb+ncoset(lb_min-1)+1,nb+ncoset(lb_max)
DO i=na+ncoset(la_min-1)+1,na+ncoset(la_max)
dabdx(i,j) = 0.0_wp
dabdy(i,j) = 0.0_wp
dabdz(i,j) = 0.0_wp
END DO
END DO
nb = nb + ncoset(lb_max)
CYCLE
END IF
! *** [da/dRi|O|b] = 2*zeta*[a+1i|O|b] - Ni(a)[a-1i|O|b] ***
DO la=MAX(0,la_min),la_max
IF (la == 0) THEN
coa = na + 1
coapx = nda + 2
coapy = nda + 3
coapz = nda + 4
DO lb=lb_min,lb_max
DO bx=0,lb
DO by=0,lb-bx
bz = lb - bx - by
cob = nb + coset(bx,by,bz)
dabdx(coa,cob) = fa*ab(coapx,cob)
dabdy(coa,cob) = fa*ab(coapy,cob)
dabdz(coa,cob) = fa*ab(coapz,cob)
END DO
END DO
END DO
ELSE
DO ax=0,la
DO ay=0,la-ax
az = la - ax - ay
coa = na + coset(ax,ay,az)
coamx = nda + coset(MAX(0,ax-1),ay,az)
coamy = nda + coset(ax,MAX(0,ay-1),az)
coamz = nda + coset(ax,ay,MAX(0,az-1))
coapx = nda + coset(ax+1,ay,az)
coapy = nda + coset(ax,ay+1,az)
coapz = nda + coset(ax,ay,az+1)
fx = REAL(ax,wp)
fy = REAL(ay,wp)
fz = REAL(az,wp)
DO lb=lb_min,lb_max
DO bx=0,lb
DO by=0,lb-bx
bz = lb - bx - by
cob = nb + coset(bx,by,bz)
dabdx(coa,cob) = fa*ab(coapx,cob) - fx*ab(coamx,cob)
dabdy(coa,cob) = fa*ab(coapy,cob) - fy*ab(coamy,cob)
dabdz(coa,cob) = fa*ab(coapz,cob) - fz*ab(coamz,cob)
END DO
END DO
END DO
END DO
END DO
END IF
END DO
nb = nb + ncoset(lb_max)
END DO
na = na + ncoset(la_max)
nda = nda + ncoset(la_max+1)
END DO
END SUBROUTINE dabdr
! *****************************************************************************
END MODULE ai_derivatives

250
src/empirical_parameters.F Normal file
View file

@ -0,0 +1,250 @@
!-----------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!!***** cp2k/empirical_parameters [1.0] *
!!
!! NAME
!! empirical_parameters
!!
!! FUNCTION
!! defines the empirical parameters for the polarization code
!! (first order hohemberg-kohn kernel and hardness kernel parameters)
!! reads them from the set file
!!
!! AUTHOR
!! gloria
!!
!! MODIFICATION HISTORY
!!
!! SOURCE
!******************************************************************************
MODULE empirical_parameters
USE atomic_kinds, ONLY : kind_info_type
USE global_types, ONLY : global_environment_type
USE kinds, ONLY : dbl
USE message_passing, ONLY : mp_bcast
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, str_search
USE termination, ONLY : stop_memory, stop_program
IMPLICIT NONE
PRIVATE
PUBLIC :: read_empirical_parameters, empirical_parameter_type
TYPE empirical_parameter_type
CHARACTER (len = 4) :: aname
INTEGER :: nset
INTEGER, DIMENSION (:), pointer :: l
REAL (dbl), DIMENSION (:), POINTER :: hk_param
REAL (dbl), DIMENSION (:), POINTER :: hardness_param
END TYPE empirical_parameter_type
!!*****
!******************************************************************************
CONTAINS
!******************************************************************************
SUBROUTINE read_empirical_parameters ( ki, set_fn, atom_names, empparm, globenv )
IMPLICIT NONE
TYPE (kind_info_type), intent(in), dimension (:) :: ki
CHARACTER ( LEN = * ), INTENT ( IN ) :: set_fn
CHARACTER ( len = * ), DIMENSION ( : ), INTENT ( IN ) :: atom_names
TYPE (empirical_parameter_type), POINTER, DIMENSION (:) :: empparm
TYPE ( global_environment_type ), INTENT ( IN ) :: globenv
! LOCALS
INTEGER :: ierror, ilen, iw, source, group, icount, nmol_type, ios, i
INTEGER :: nkind, nset, iset, iat, natom_types
CHARACTER ( LEN = 20 ) :: string, string2
CHARACTER ( LEN = 6 ) :: label, aname, merda, shit
CHARACTER ( LEN = 6 ), POINTER, DIMENSION (:) :: element_symbol
TYPE (empirical_parameter_type), POINTER, DIMENSION (:) :: read_ep
iw = globenv % scr
icount = 0
nkind = 0
ALLOCATE (read_ep(size(ki)), stat = ios)
IF (ios /= 0 ) CALL stop_memory &
( 'read_empirical_parameters', 'read_ep', nkind )
ALLOCATE (empparm(size(ki)), stat = ios)
IF (ios /= 0 ) CALL stop_memory &
( 'read_empirical_parameters', 'empparm', nkind )
ALLOCATE (element_symbol(size(ki)), stat = ios)
IF (ios /= 0 ) CALL stop_memory &
( 'read_emporical_parameters', 'elementsymbol', nkind )
!..parse the input section
label = '&EMPIRICAL'
CALL parser_init ( set_fn, label, ierror, globenv )
IF ( ierror /= 0 ) THEN
IF ( globenv % ionode ) THEN
WRITE ( iw, '( A )' ) ' No input section &EMPIRICAL found on file '
WRITE ( iw, '( T2, A )' ) set_fn
END IF
CALL stop_parser ( 'read_empirical_parameters', '&EMPIRICAL' )
ELSE
CALL read_line
DO WHILE (test_next()/='X')
ilen = 8
CALL cfield(string,ilen)
natom_types = SIZE ( atom_names )
iat = str_search(atom_names,natom_types,string)
nkind = nkind + 1
DO
CALL read_line
ilen = 8
CALL cfield(string2,ilen)
CALL uppercase ( string2 )
SELECT CASE (string2)
CASE DEFAULT
CALL p_error()
CALL stop_parser ( 'read_empirical_parameter','unknown option')
CASE ( 'NSET')
nset = get_int()
read_ep (iat) % nset = nset
IF ( read_ep (iat) % nset > 0) THEN
ALLOCATE (read_ep(iat)%l(nset), STAT = ios)
ALLOCATE (read_ep(iat)%hk_param(nset), stat = ios)
ALLOCATE (read_ep(iat)%hardness_param(nset), stat = ios)
IF (ios /= 0 ) CALL stop_memory &
( 'read_emporical_parameters', 'read_ep%l', nset )
element_symbol (iat) = string
END IF
CASE ( 'EMPAR')
i = get_int ()
read_ep(iat) % l (i) = get_int()
read_ep(iat) % hk_param (i) = get_real()
read_ep(iat) % hardness_param (i) = get_real()
CASE ( 'END')
ilen = 0
CALL cfield(string2,ilen)
EXIT
END SELECT
END DO
CALL read_line
END DO
END IF
CALL parser_end
! map the empirical parameter to the corresponding structure
! and checks that the mapping is the same as in kind_info
IF (nkind /= size(ki)) CALL stop_program ('read_empirical_parameters', &
'inconsistent number of polarizable atomic kinds')
DO i = 1, nkind
aname = ki(i)% element_symbol
iat = str_search(element_symbol,nkind,aname)
empparm (i) % aname = aname
nset = read_ep(iat)% nset
empparm (i) % nset = read_ep(iat)% nset
IF (empparm(i)% nset /= ki(i) % orb_basis_set % nset) CALL stop_program &
('read_empirical_parameters', 'number of polarizable sets not matching')
ALLOCATE ( empparm(i) % l (nset), STAT = ios )
ALLOCATE ( empparm(i) % hk_param (nset), stat = ios )
ALLOCATE ( empparm(i) % hardness_param (nset), stat = ios )
IF (ios /= 0 ) CALL stop_memory &
( 'read_molecule_section', 'empparm%l', nset )
DO iset = 1, nset
empparm(i) % l (iset) = read_ep ( iat ) % l (iset)
empparm(i) % hk_param (iset) = read_ep ( iat ) % hk_param (iset)
empparm(i) % hardness_param (iset) = read_ep ( iat ) % hardness_param (iset)
END DO
END DO
! ..write some information to output
IF (globenv%ionode) THEN
IF (globenv%print_level>=0) THEN
DO i = 1, nkind
WRITE ( iw, '( A,T61,A )' ) ' EMPIRICAL_PARAMETERS| atom type ', &
ADJUSTR ( empparm(i) % aname )
WRITE ( iw, '( A,T71,I10 )' ) ' EMPIRICAL_PARAMETERS| Number of sets ', &
empparm(i) % nset
DO iset = 1, empparm (i) % nset
WRITE ( iw, '( A,T71,I10 )' ) ' EMPIRICAL| angular momentum number ', &
empparm(i) % l (iset)
WRITE ( iw, '( A,T71,F12.6 )' ) ' EMPIRICAL| Hohemberg-Kohn parameter ', &
empparm(i) % hk_param (iset)
WRITE ( iw, '( A,T71,F12.6 )' ) ' EMPIRICAL| Hardness parameter ', &
empparm(i) % hardness_param (iset)
END DO
WRITE ( iw,'( )' )
END DO
END IF
END IF
! deallocates
DO i = 1, nkind
IF (ASSOCIATED (read_ep (i) %l)) DEALLOCATE (read_ep (i) %l, stat = ios)
IF (ios /= 0 ) CALL stop_memory &
( 'read_molecule_section', ' deallocate read_ep%l' )
IF (ASSOCIATED(read_ep (i) %hk_param)) DEALLOCATE (read_ep (i) %hk_param, stat = ios)
IF (ios /= 0 ) CALL stop_memory &
( 'read_molecule_section', ' deallocate read_ep%hk_param' )
IF (ASSOCIATED(read_ep (i) %hardness_param )) &
DEALLOCATE (read_ep (i) % hardness_param, stat = ios)
IF (ios /= 0 ) CALL stop_memory &
( 'read_molecule_section', ' deallocate read_ep%l' )
END DO
IF (ASSOCIATED(read_ep)) DEALLOCATE (read_ep, stat = ios)
IF (ios /= 0 ) CALL stop_memory &
( 'read_molecule_section', ' deallocate read_ep' )
END SUBROUTINE read_empirical_parameters
!!*****
!******************************************************************************
END MODULE empirical_parameters
!!*****
!******************************************************************************

756
src/pol_overlap.F Normal file
View file

@ -0,0 +1,756 @@
!-----------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 1999 MPI fuer Festkoerperforschung, Stuttgart !
!-----------------------------------------------------------------------------!
!!****** cp2k/pol_overlap [1.0] *
!!
!! NAME
!! pol_overlap
!!
!! FUNCTION
!! Calculation of the second order contribution of the Hohenbeg-Kohn
!! functional that depends on the overlaps between basis functions
!!
!! E (overl) = 1/2 sum_ij c_i*c_j (f_i + f_j) <phi_i|phi_j>
!!
!! AUTHOR
!! gloria (06 Jun 2001)
!!
!! MODIFICATION HISTORY
!! none
!!
!! SOURCE
!******************************************************************************
MODULE pol_overlap
! *****************************************************************************
USE kinds, ONLY: dbl
USE ai_derivatives, ONLY : dabdr
USE ai_overlap, ONLY : overlap
USE ao_types, ONLY : ao_type
USE atomic_kinds, ONLY: kind_info_type
USE empirical_parameters, ONLY : empirical_parameter_type
USE fist_nonbond_force, ONLY : find_image
USE global_types, ONLY: global_environment_type
USE mathlib, ONLY: symmetrize_matrix
USE matrix_types, ONLY: add_block_node,&
allocate_matrix,&
deallocate_matrix,&
deallocate_matrix_row,&
first_block_node,&
get_block_node,&
next_block_node,&
real_block_node_type,&
real_matrix_set_type
USE md, ONLY : thermodynamic_type
USE memory_utilities, ONLY: reallocate
USE method_specifications, ONLY: maxder
USE molecule_types, ONLY : particle_node_type, &
linklist_images, linklist_neighbor
USE orbital_pointers, ONLY: nco,ncoset
USE particle_types, ONLY : particle_type
USE simulation_cell, ONLY : cell_type, get_cell_param
USE termination, ONLY: stop_memory, stop_program
USE timings, ONLY: timeset,timestop
IMPLICIT NONE
PRIVATE
PUBLIC :: force_overlap
TYPE atom_basis_info_type
INTEGER :: ipart, ikind
INTEGER :: first_cgf, last_cgf
END TYPE atom_basis_info_type
! *****************************************************************************
CONTAINS
! *****************************************************************************
SUBROUTINE force_overlap(ao, kind_info, part, pnode, box, thermo, empparm, &
fo_coef, fo_part, globenv)
! ***************************************************************************
TYPE ( ao_type ), INTENT(INOUT) :: ao
TYPE ( kind_info_type ), INTENT (IN), DIMENSION (:) :: kind_info
TYPE ( particle_type ), INTENT (IN), dimension (:) :: part
TYPE ( particle_node_type ), INTENT (INOUT), dimension (:) :: pnode
TYPE ( cell_type ), INTENT ( INOUT ) :: box
TYPE ( thermodynamic_type ), INTENT ( INOUT ) :: thermo
TYPE ( empirical_parameter_type ), DIMENSION ( : ), INTENT ( IN ) :: empparm
REAL ( dbl ), DIMENSION (:), INTENT(INOUT) :: fo_coef
REAL ( dbl ), DIMENSION (:,:), INTENT(INOUT), OPTIONAL :: fo_part
TYPE(global_environment_type), INTENT(IN) :: globenv
! locals
INTEGER :: maxco,maxcgf,maxdim
REAL(dbl), DIMENSION(:,:), ALLOCATABLE :: sab,work
REAL(dbl), DIMENSION(:,:), ALLOCATABLE :: sabdx, sabdy, sabdz
REAL (dbl), DIMENSION (:,:,:), ALLOCATABLE :: atom_block
REAL(dbl) :: radius_a,radius_b, rijsq, r2max, rcut
INTEGER :: first_cgfa,first_cgfb,id,i,ii,is,im,j,ipart,jpart, &
iatom,ikind,ineighbor,ipgf,iset,ishell,istat,&
jatom,jkind,jneighbor,jpgf,jset,jshell,&
la_max,la_min,last_cgfa,last_cgfb,lb_max,lb_min,&
ncgfa,ncgfb,ncoa,ncob,npgfa,npgfb,nseta,nsetb,nshella,nshellb,&
ncgf,natom,nat,npart,nkind
REAL (dbl), DIMENSION (:,:), POINTER :: s_block
REAL (dbl), DIMENSION (:), POINTER :: zeta, zetb, rpgfa, rpgfb
REAL (dbl) :: energy
TYPE (real_matrix_set_type) :: smat
REAL ( dbl ), DIMENSION (3) :: ri, rij, s, perd, &
vec, quotient, cell_lengths
REAL ( dbl ), DIMENSION (3,3) :: hmat, h_inv
TYPE (linklist_neighbor), POINTER :: current_neighbor
TYPE (linklist_images), POINTER :: current_image
TYPE (atom_basis_info_type), DIMENSION (:), POINTER :: ainfo
INTEGER, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: n_images
LOGICAL :: first_time
LOGICAL :: forces
! ---------------------------------------------------------------------------
! initialize atomic info
ncgf = 0
natom = 0
nat = 0
npart = size(pnode)
nkind = size(kind_info)
hmat = box % hmat
h_inv = box % h_inv
perd = box % perd
DO ikind = 1, nkind
nat = nat + size (kind_info(ikind) % atom_list)
END DO
NULLIFY (ainfo)
IF (.NOT. ASSOCIATED (ainfo)) then
ALLOCATE (ainfo(nat), stat=istat)
IF (istat /= 0) THEN
CALL stop_memory("pol_overlap","ainfo",nat)
END IF
END IF
DO i = 1,npart
ikind = pnode (i) % p % kind
if (ikind == 0) CYCLE
natom = natom + 1
pnode ( i ) % p % abase = natom
ainfo ( natom ) % ipart = i
ainfo ( natom ) % ikind = ikind
pnode ( i ) % p % first_cgf = ncgf + 1
ainfo ( natom ) % first_cgf = pnode ( i ) % p % first_cgf
ncgf = ncgf + kind_info(ikind)%orb_basis_set%ncgf
pnode ( i ) % p % last_cgf = ncgf
ainfo ( natom ) % last_cgf = pnode ( i ) % p % last_cgf
END DO
if (nat /= natom) call stop_program("pol_overlap","inconsistent number of atoms with basis")
! Allocate work storage
maxco = 0
maxcgf = 0
IF (maxder==1) THEN
maxdim = 4
forces = .TRUE.
ELSE IF (maxder==0) THEN
maxdim = 1
forces = .FALSE.
END If
DO ikind=1,nkind
DO iset=1,kind_info(ikind)%orb_basis_set%nset
maxcgf = MAX(maxcgf,kind_info(ikind)%orb_basis_set%ncgf)
npgfa = kind_info(ikind)%orb_basis_set%npgf(iset)
la_max = kind_info(ikind)%orb_basis_set%lmax(iset) + maxder
maxco = MAX(maxco,npgfa*ncoset(la_max))
END DO
END DO
ALLOCATE (sab(maxco,maxco),STAT=istat)
IF (istat /= 0) THEN
CALL stop_memory("pol_overlap","sab",maxco*maxco)
END IF
sab(:,:) = 0.0_dbl
ALLOCATE (sabdx(maxco,maxco),STAT=istat)
IF (istat /= 0) THEN
CALL stop_memory("pol_overlap","sabdx",maxco*maxco)
END IF
sabdx(:,:) = 0.0_dbl
ALLOCATE (sabdy(maxco,maxco),STAT=istat)
IF (istat /= 0) THEN
CALL stop_memory("pol_overlap","sabdy",maxco*maxco)
END IF
sabdy(:,:) = 0.0_dbl
ALLOCATE (sabdz(maxco,maxco),STAT=istat)
IF (istat /= 0) THEN
CALL stop_memory("pol_overlap","sabdz",maxco*maxco)
END IF
sabdz(:,:) = 0.0_dbl
ALLOCATE (work(maxco,maxcgf),STAT=istat)
IF (istat /= 0) THEN
CALL stop_memory("pol_overlap","work",maxco*maxcgf)
END IF
work(:,:) = 0.0_dbl
ALLOCATE (atom_block(maxco,maxcgf,maxdim),STAT=istat)
IF (istat /= 0) THEN
CALL stop_memory("pol_overlap","atom_blok",maxco*maxcgf*maxdim)
END IF
atom_block(:,:,:) = 0.0_dbl
NULLIFY (smat%matrix)
CALL build_block_overlap_matrix(pnode,part,ainfo,smat,natom,ncgf)
! initialize forces
fo_coef = 0._dbl
IF ( PRESENT ( fo_part ) ) fo_part = 0._dbl
! intialize images
first_time = .NOT.ALLOCATED ( n_images )
IF ( first_time ) THEN
IF ( .NOT. ALLOCATED ( n_images ) ) &
ALLOCATE ( n_images ( nkind, nkind, 3 ), STAT = istat )
IF ( istat /= 0 ) CALL stop_memory ( 'force_control', &
'n_images', nkind ** 2 * 3 )
CALL get_cell_param ( box, cell_lengths )
DO ikind = 1, nkind
DO jkind = 1, nkind
rcut = kind_info( ikind ) % orb_basis_set % kind_radius + &
kind_info( jkind) % orb_basis_set % kind_radius
quotient ( : ) = rcut / cell_lengths ( : )
DO id=1,3
IF ( quotient ( id ) <= 0.5_dbl ) THEN
n_images ( ikind, jkind, id ) = 0
ELSE
n_images ( ikind, jkind, id ) &
= CEILING ( rcut / cell_lengths ( id ) )
END IF
END DO
END DO
END DO
END IF
! Loop over all atoms
DO i = 1, size(pnode)
! cycle if the atoms has no basis function!
ipart = pnode(i) %p%iatom
ikind = part (ipart) % kind
IF (ikind <= 0) CYCLE
nseta = kind_info(ikind)%orb_basis_set%nset
iatom = pnode (i) % p % abase
ri = pnode (i) % p % r
! Loop over all neighbor atoms of the current atom "iatom"
current_neighbor => pnode ( i ) % sl
DO j = 1, pnode (i) % nsneighbor
jpart = current_neighbor % index
jatom = part (jpart) % abase
rij = current_neighbor % p % r - ri
jkind = current_neighbor % p % kind
nsetb = kind_info(jkind)%orb_basis_set%nset
DO iset=1,nseta
la_max = kind_info(ikind)%orb_basis_set%lmax(iset)
la_min = kind_info(ikind)%orb_basis_set%lmin(iset)
npgfa = kind_info(ikind)%orb_basis_set%npgf(iset)
nshella = kind_info(ikind)%orb_basis_set%nshell(iset)
rpgfa => kind_info(ikind)%orb_basis_set%pgf_radius(1:npgfa,iset)
zeta => kind_info(ikind)%orb_basis_set%zet(1:npgfa,iset)
first_cgfa = kind_info(ikind)%orb_basis_set%first_cgf(1,iset)
last_cgfa = kind_info(ikind)%orb_basis_set%last_cgf(nshella,iset)
ncgfa = last_cgfa - first_cgfa + 1
ncoa = npgfa*ncoset(la_max)
IF (iatom == jatom) THEN
is = iset
ELSE
is = 1
END IF
radius_a = kind_info(ikind)%orb_basis_set%set_radius(iset)
DO jset=is,nsetb
radius_b = kind_info(jkind)%orb_basis_set%set_radius(jset)
s(1) = h_inv(1,1)*rij(1) + h_inv(1,2)*rij(2) + h_inv(1,3)*rij(3)
s(2) = h_inv(2,1)*rij(1) + h_inv(2,2)*rij(2) + h_inv(2,3)*rij(3)
s(3) = h_inv(3,1)*rij(1) + h_inv(3,2)*rij(2) + h_inv(3,3)*rij(3)
IF ( MAXVAL ( n_images ( ikind, jkind, : ) ) == 0 ) THEN
IF ( ABS ( s ( 1 ) ) > 0.5_dbl ) THEN
s(1) = s(1) - perd(1) * INT(s(1)+SIGN(0.5_dbl,s(1)))
END IF
IF ( ABS ( s ( 2 ) ) > 0.5_dbl ) THEN
s(2) = s(2) - perd(2) * INT(s(2)+SIGN(0.5_dbl,s(2)))
END IF
IF ( ABS ( s ( 3 ) ) > 0.5_dbl ) THEN
s(3) = s(3) - perd(3) * INT(s(3)+SIGN(0.5_dbl,s(3)))
END IF
rij(1) = hmat(1,1)*s(1) + hmat(1,2)*s(2) + hmat(1,3)*s(3)
rij(2) = hmat(2,1)*s(1) + hmat(2,2)*s(2) + hmat(2,3)*s(3)
rij(3) = hmat(3,1)*s(1) + hmat(3,2)*s(2) + hmat(3,3)*s(3)
END IF
rijsq = rij ( 1 ) ** 2 + rij ( 2 ) ** 2 + rij ( 3 ) ** 2
r2max = (radius_a + radius_b)**2
IF (rijsq > r2max ) CYCLE
CALL get_block_node(matrix=smat%matrix,&
block_row=iatom,&
block_col=jatom,&
block=s_block)
IF (ASSOCIATED(s_block)) THEN
atom_block(:,:,:) = 0.0_dbl
lb_max = kind_info(jkind)%orb_basis_set%lmax(jset)
lb_min = kind_info(jkind)%orb_basis_set%lmin(jset)
npgfb = kind_info(jkind)%orb_basis_set%npgf(jset)
nshellb = kind_info(jkind)%orb_basis_set%nshell(jset)
rpgfb => kind_info(jkind)%orb_basis_set%pgf_radius(1:npgfb,jset)
zetb => kind_info(jkind)%orb_basis_set%zet(1:npgfb,jset)
first_cgfb = kind_info(jkind)%orb_basis_set%first_cgf(1,jset)
last_cgfb = kind_info(jkind)%orb_basis_set%last_cgf(nshellb,jset)
ncgfb = last_cgfb - first_cgfb + 1
ncob = npgfb*ncoset(lb_max)
! Calculate the overlap integrals and the forces on the central image
CALL get_overlap (la_max,zeta,rpgfa,la_min,lb_max,zetb,rpgfb,lb_min, &
rij,rijsq,ncoa,ncob,ncgfa,ncgfb,sab,sabdx,sabdy, &
sabdz,work,s_block(first_cgfa:last_cgfa,first_cgfb:last_cgfb), &
kind_info(ikind)%orb_basis_set%cphi(1:ncoa,first_cgfa:last_cgfa), &
kind_info(jkind)%orb_basis_set%cphi(1:ncob,first_cgfb:last_cgfb), &
atom_block(first_cgfa:last_cgfa,first_cgfb:last_cgfb,1:maxdim), forces)
! summing over lattice translations of neighbors
current_image => current_neighbor % image
DO im = 1, current_neighbor % nimages
CALL find_image ( s, perd, current_image % vec, hmat, rijsq, rij )
IF ( rijsq <= r2max ) THEN
CALL get_overlap (la_max,zeta,rpgfa,la_min,lb_max,zetb,rpgfb,lb_min, &
rij,rijsq,ncoa,ncob,ncgfa,ncgfb,sab,sabdx,sabdy, &
sabdz,work,s_block(first_cgfa:last_cgfa,first_cgfb:last_cgfb), &
kind_info(ikind)%orb_basis_set%cphi(1:ncoa,first_cgfa:last_cgfa), &
kind_info(jkind)%orb_basis_set%cphi(1:ncob,first_cgfb:last_cgfb), &
atom_block(first_cgfa:last_cgfa,first_cgfb:last_cgfb,1:maxdim), forces)
END IF
current_image => current_image % next
END DO
END IF
! here the contribution to the forces on atoms and coeffs given
! from each atom_block is calculated
CALL get_forces_on_coefs(ao,fo_coef,energy,empparm,atom_block, &
part,ipart,jpart,iset,jset)
IF (forces.AND.(iatom/=jatom)) CALL get_forces_on_atoms (ao, fo_part, &
empparm, atom_block, part, ipart, jpart, iset, jset)
END DO
END DO
! Symmetrize the diagonal blocks ***
IF (iatom == jatom) THEN
CALL symmetrize_matrix(s_block,"upper_to_lower")
END IF
current_neighbor => current_neighbor % next
END DO
END DO
IF (ALLOCATED(sab)) THEN
DEALLOCATE (sab,STAT=istat)
IF (istat /= 0) THEN
CALL stop_memory("pol_overlap","sab")
END IF
END IF
IF (ALLOCATED(sabdx)) THEN
DEALLOCATE (sabdx,STAT=istat)
IF (istat /= 0) THEN
CALL stop_memory("pol_overlap","sabdx")
END IF
END IF
IF (ALLOCATED(sabdy)) THEN
DEALLOCATE (sabdy,STAT=istat)
IF (istat /= 0) THEN
CALL stop_memory("pol_overlap","sabdy")
END IF
END IF
IF (ALLOCATED(sabdz)) THEN
DEALLOCATE (sabdz,STAT=istat)
IF (istat /= 0) THEN
CALL stop_memory("pol_overlap","sabdz")
END IF
END IF
IF (ALLOCATED(work)) THEN
DEALLOCATE (work,STAT=istat)
IF (istat /= 0) THEN
CALL stop_memory("pol_overlap","work")
END IF
END IF
IF (ALLOCATED(atom_block)) THEN
DEALLOCATE (atom_block,STAT=istat)
IF (istat /= 0) THEN
CALL stop_memory("pol_overlap","atom_block")
END IF
END IF
IF (ASSOCIATED(ainfo)) THEN
DEALLOCATE (ainfo,STAT=istat)
IF (istat /= 0) THEN
CALL stop_memory("pol_overlap","ainfo")
END IF
END IF
CALL deallocate_matrix(smat%matrix)
END SUBROUTINE force_overlap
!---------------------------------------------------------------------------------
SUBROUTINE get_overlap(la_max,zeta,rpgfa,la_min,lb_max,zetb,rpgfb,lb_min, &
rij,rijsq,ncoa,ncob,ncgfa,ncgfb,sab,sabdx,sabdy, &
sabdz,work,s_block,cphia,cphib,atom_block,forces)
!---------------------------------------------------------------------------------
IMPLICIT NONE
INTEGER, INTENT (IN) :: la_max, la_min, lb_max, lb_min
REAL (dbl), DIMENSION (:), POINTER :: rpgfa,rpgfb,zeta,zetb
INTEGER, INTENT (IN) :: ncoa, ncob, ncgfa, ncgfb
REAL (dbl), INTENT (IN), DIMENSION (:) :: rij
REAL (dbl), intent (IN) :: rijsq
REAL(dbl), DIMENSION (:,:), intent(in) :: cphia, cphib
REAL(dbl), DIMENSION (:,:), intent(inout) :: sab, sabdx, sabdy, sabdz
REAL(dbl), DIMENSION (:,:), intent(inout) :: s_block
REAL(dbl), DIMENSION (:,:), intent(inout) :: work
REAL(dbl), DIMENSION (:,:,:), intent (inout) :: atom_block
INTEGER :: lda, ldb, ldc, ldd
LOGICAL, INTENT (IN) :: forces
! locals
REAL (dbl) :: dab
! ---------------------------------------------------------------------------
dab = sqrt(rijsq)
! calculate the primitive overlap integral
CALL overlap(la_max,zeta,rpgfa,la_min,&
lb_max,zetb,rpgfb,lb_min,&
rij,rijsq,sab)
! Contraction step (overlap matrix)
lda = size(sab,1)
ldb = size(cphib,1)
ldc = size(work,1)
CALL dgemm("N","N",ncoa,ldd,ncob,1.0_dbl,sab(1,1),lda,&
cphib(1,1), &
ldb,0.0_dbl,&
work(1,1),ldc)
lda = size(cphia,1)
ldb = size(work,1)
ldc = size(s_block,1)
CALL dgemm("T","N",ncgfa,ncgfb,ncoa,1.0_dbl,&
cphia,&
lda,&
work,ldb,1.0_dbl,&
s_block,ldc)
! get the total overlap integral
atom_block(:,:,1) = atom_block(:,:,1) + s_block(:,:)
! if requested, calculate the derivatives
if (forces) then
CALL dabdr(la_max,zeta,rpgfa,la_min,&
lb_max,zetb,rpgfb,lb_min,&
dab,sab,sabdx,sabdy,sabdz)
! contraction step for the x component
CALL dgemm("N","N",ncoa,ncgfb,ncob,1.0_dbl,sabdx(1,1),SIZE(sabdx,1),&
cphib, SIZE(cphib,1),0.0_dbl,work(1,1),SIZE(work,1))
CALL dgemm("T","N",ncgfa,ncgfb,ncoa,1.0_dbl,cphia,SIZE(cphia,1),&
work(1,1),SIZE(work,1),1.0_dbl,s_block,SIZE(s_block,1))
! get the total x derivative
atom_block(:,:,2) = atom_block(:,:,2) + s_block(:,:)
! contraction step for the y component
CALL dgemm("N","N",ncoa,ncgfb,ncob,1.0_dbl,sabdy(1,1),SIZE(sabdy,1),&
cphib, SIZE(cphib,1),0.0_dbl,work(1,1),SIZE(work,1))
CALL dgemm("T","N",ncgfa,ncgfb,ncoa,1.0_dbl,cphia,SIZE(cphia,1),&
work(1,1),SIZE(work,1),1.0_dbl,s_block,SIZE(s_block,1))
! get the total y derivative
atom_block(:,:,3) = atom_block(:,:,3) + s_block(:,:)
! contraction step for the z component
CALL dgemm("N","N",ncoa,ncgfb,ncob,1.0_dbl,sabdz(1,1),SIZE(sabdz,1),&
cphib, SIZE(cphib,1),0.0_dbl,work(1,1),SIZE(work,1))
CALL dgemm("T","N",ncgfa,ncgfb,ncoa,1.0_dbl,cphia,SIZE(cphia,1),&
work(1,1),SIZE(work,1),1.0_dbl,s_block,SIZE(s_block,1))
! get the total z derivative
atom_block(:,:,4) = atom_block(:,:,4) + s_block(:,:)
END IF
END SUBROUTINE get_overlap
!--------------------------------------------------------------------------------------
SUBROUTINE build_block_overlap_matrix(pnode,part,ainfo,smat,natom,ncgf)
!--------------------------------------------------------------------------------------
IMPLICIT NONE
TYPE (particle_node_type), dimension (:), intent (in) :: pnode
TYPE (particle_type), dimension (:), intent (in) :: part
TYPE (atom_basis_info_type), DIMENSION (:) , INTENT (IN) :: ainfo
TYPE (real_matrix_set_type), intent (out) :: smat
INTEGER, intent (in) :: natom
INTEGER, intent (in) :: ncgf
! locals
INTEGER :: i, ikind, ii, ipart, iatom, j, jpart, jatom
TYPE (linklist_neighbor), POINTER :: current_neighbor
CALL allocate_matrix(matrix=smat%matrix,&
nblock_row=natom,&
nblock_col=natom,&
nrow=ncgf,&
ncol=ncgf,&
first_row=ainfo(:)%first_cgf,&
last_row=ainfo(:)%last_cgf,&
first_col=ainfo(:)%first_cgf,&
last_col=ainfo(:)%last_cgf,&
matrix_name="OVERLAP MATRIX",&
matrix_symmetry="symmetric")
DO i = 1, size(pnode)
ipart = pnode(i) %p%iatom
ikind = part (ipart) % kind
! cycle if the particle has no polarization basis functions
IF (ikind <= 0) CYCLE
iatom = pnode (i) % p % abase
! Loop over all neighbor atoms of the current atom "iatom"
! (the neighbor list is constructed with the black-white scheme
! and the blocks are allocated accordingly)
current_neighbor => pnode ( i ) % sl
DO j = 1, pnode (i) % nsneighbor
jpart = current_neighbor % index
jatom = part (jpart) % abase
CALL add_block_node ( matrix=smat%matrix, block_row=iatom, block_col=jatom )
END DO
current_neighbor => current_neighbor % next
END DO
END SUBROUTINE build_block_overlap_matrix
!--------------------------------------------------------------------------------------
SUBROUTINE get_forces_on_coefs(ao,fo_coef,energy,empparm,atom_block, &
part,ipart,jpart,iset,jset)
!--------------------------------------------------------------------------------------
IMPLICIT NONE
TYPE (ao_type), intent (inout) :: ao
REAL (dbl), DIMENSION (:), INTENT(INOUT) :: fo_coef
REAL (dbl), INTENT (INOUT) :: energy
TYPE ( empirical_parameter_type ), DIMENSION ( : ), INTENT ( IN ) :: empparm
REAL (dbl), INTENT (IN) , DIMENSION (:,:,:) :: atom_block
TYPE (particle_type), dimension (:), intent (in) :: part
INTEGER, intent (in) :: ipart,jpart
INTEGER, intent (in) :: iset,jset
! locals
INTEGER :: ikind, jkind, first_cgfa, last_cgfa, first_cgfb, last_cgfb
INTEGER :: nshellb, nshella, ncgfa, ncgfb, ic, jc, icoef, jcoef
REAL ( dbl ) :: hpi, hpj
ikind = part (ipart) % kind
jkind = part (jpart) % kind
first_cgfa = ao % kind_info (ikind) % orb_basis_set % first_cgf (1,iset)
nshella = ao % kind_info (ikind) % orb_basis_set % nshell (iset)
last_cgfa = ao % kind_info (ikind) % orb_basis_set % last_cgf (nshella,iset)
first_cgfb = ao % kind_info (jkind) % orb_basis_set % first_cgf (1,jset)
nshellb = ao % kind_info (jkind) % orb_basis_set % nshell (jset)
last_cgfb = ao % kind_info (jkind) % orb_basis_set % last_cgf (nshellb,jset)
hpi = empparm (ikind) % hardness_param (iset)
hpj = empparm (jkind) % hardness_param (jset)
DO ic = first_cgfa, last_cgfa
icoef = part (ipart) % coef_list (ic)
DO jc = first_cgfb, last_cgfb
jcoef = part (jpart) % coef_list (jc)
fo_coef (icoef) = fo_coef (icoef) + (hpi + hpj) * ao % cr (jcoef) * &
atom_block (ic,jc,1)
fo_coef (jcoef) = fo_coef (jcoef) + (hpi + hpj) * ao % cr (icoef) * &
atom_block (ic,jc,1)
energy = energy + 0.5 * (hpi + hpj) * ao % cr (icoef) * &
ao % cr (jcoef) * atom_block (ic, jc, 1)
END DO
END DO
END SUBROUTINE get_forces_on_coefs
!--------------------------------------------------------------------------------------
SUBROUTINE get_forces_on_atoms(ao,fo_part,empparm,atom_block, &
part,ipart,jpart,iset,jset)
!--------------------------------------------------------------------------------------
IMPLICIT NONE
TYPE (ao_type), intent (inout) :: ao
REAL (dbl), DIMENSION (:,:), INTENT(INOUT) :: fo_part
TYPE ( empirical_parameter_type ), DIMENSION ( : ), INTENT ( IN ) :: empparm
REAL (dbl), INTENT (IN) , DIMENSION (:,:,:) :: atom_block
TYPE (particle_type), dimension (:), intent (in) :: part
INTEGER, intent (in) :: ipart,jpart
INTEGER, intent (in) :: iset,jset
! locals
INTEGER :: ikind, jkind, first_cgfa, last_cgfa, first_cgfb, last_cgfb
INTEGER :: nshellb, nshella, ncgfa, ncgfb, ic, jc, icoef, jcoef
REAL ( dbl ) :: hpi, hpj
ikind = part (ipart) % kind
jkind = part (jpart) % kind
first_cgfa = ao % kind_info (ikind) % orb_basis_set % first_cgf (1,iset)
nshella = ao % kind_info (ikind) % orb_basis_set % nshell (iset)
last_cgfa = ao % kind_info (ikind) % orb_basis_set % last_cgf (nshella,iset)
first_cgfb = ao % kind_info (jkind) % orb_basis_set % first_cgf (1,jset)
nshellb = ao % kind_info (jkind) % orb_basis_set % nshell (jset)
last_cgfb = ao % kind_info (jkind) % orb_basis_set % last_cgf (nshellb,jset)
hpi = empparm (ikind) % hardness_param (iset)
hpj = empparm (jkind) % hardness_param (jset)
DO ic = first_cgfa, last_cgfa
icoef = part (ipart) % coef_list (ic)
DO jc = first_cgfb, last_cgfb
jcoef = part (jpart) % coef_list (jc)
fo_part (1,ipart) = fo_part (1,ipart) + (hpi + hpj) * ao % cr (jcoef) * &
atom_block (ic,jc,2)
fo_part (1,jpart) = fo_part (1,jpart) - (hpi + hpj) * ao % cr (icoef) * &
atom_block (ic,jc,2)
fo_part (2,ipart) = fo_part (2,ipart) + (hpi + hpj) * ao % cr (jcoef) * &
atom_block (ic,jc,3)
fo_part (2,jpart) = fo_part (2,jpart) - (hpi + hpj) * ao % cr (icoef) * &
atom_block (ic,jc,3)
fo_part (3,ipart) = fo_part (3,ipart) + (hpi + hpj) * ao % cr (jcoef) * &
atom_block (ic,jc,4)
fo_part (3,jpart) = fo_part (3,jpart) - (hpi + hpj) * ao % cr (icoef) * &
atom_block (ic,jc,4)
END DO
END DO
END SUBROUTINE get_forces_on_atoms
! *****************************************************************************
END MODULE pol_overlap
! *****************************************************************************