diff --git a/src/OBJECTDEFS b/src/OBJECTDEFS index 04e9ca046e..d13c8d304b 100644 --- a/src/OBJECTDEFS +++ b/src/OBJECTDEFS @@ -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\ diff --git a/src/ai_derivatives.F b/src/ai_derivatives.F new file mode 100644 index 0000000000..1a7de19902 --- /dev/null +++ b/src/ai_derivatives.F @@ -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 + + diff --git a/src/empirical_parameters.F b/src/empirical_parameters.F new file mode 100644 index 0000000000..a70b5b66c1 --- /dev/null +++ b/src/empirical_parameters.F @@ -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 + +!!***** +!****************************************************************************** diff --git a/src/pol_overlap.F b/src/pol_overlap.F new file mode 100644 index 0000000000..03837d7972 --- /dev/null +++ b/src/pol_overlap.F @@ -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) +!! +!! 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 + +! *****************************************************************************