1module qmcmm_response
2
3   implicit none
4
5   public qmcmm_lr
6   public qmcmm_qr
7
8   private
9
10contains
11
12   subroutine qmcmm_lr(nosim,bovecs,cref,cmo,udv,dv,udvtr,   &
13                            dvtr,evecs,work,lwork)
14!
15! Purpose:
16!     Computes QM/NP/MM contribution to linear response vector.
17!
18! Input:
19!   NOSIM  - Number of first Fock/Kohn-Sham matrices
20!   BOVECs - Linear response vectors
21!   CREF   - Reference state CI coeficients (Currently not used)
22!   CMO    - Molecular orbitals
23!   UDV    - Density matrix
24!   DV     - Active density matrix
25!   UDVTR  - Triplet density matrix
26!   DVTR   - Active triplet density matrix
27!   WORK   - Temporary memory array
28!   LWORK  - Size of temporary memory array
29! Output:
30!   EVECS  - Response vectors (XY).
31!
32! Last updated: 22/03/2013 by Z. Rinkevicius.
33!
34
35      use qmcmm, only: getdim_relmat, read_relmat
36
37#include "inforb.h"
38#include "infdim.h"
39#include "wrkrsp.h"
40#include "infrsp.h"
41#include "priunit.h"
42#include "qmnpmm.h"
43
44      integer, intent(in) :: nosim
45      integer, intent(in) :: lwork
46
47      real(8) :: BOVECS(*), CMO(*), UDV(NASHDI,NASHDI)
48      real(8) :: UDVTR(N2ASHX), DVTR(*), EVECS(KZYVAR,*)
49      real(8) :: WORK(*), DV(*), CREF(*)
50      integer :: i, j, idimension, ioff
51
52      real(8), allocatable :: mqvec(:)
53      real(8), allocatable :: fvvec(:)
54      real(8), allocatable :: ucmo(:)
55      real(8), allocatable :: bov(:)
56      real(8), allocatable :: relmat(:)
57      real(8), allocatable :: rxy(:)
58      real(8), allocatable :: rxyt(:)
59
60!     np/mm contribution to response matrix is zero for triplet
61!     perturbations applied to singlet reference state
62      if ((nasht == 0) .and. (trplet)) return
63
64      allocate(ucmo(nbast*norbt))
65      ucmo = 0.0d0
66      allocate(bov(nosim*n2orbx))
67      bov = 0.0d0
68      allocate(rxy(nosim*n2orbx))
69      rxy = 0.0d0
70      if (trplet) then
71          allocate(rxyt(nosim*n2orbx))
72          rxyt = 0.0d0
73      end if
74
75      if (mqiter) then ! iterative method
76         call quit('mqiter not implemented in rspqmnp.F90')
77      else ! non iterative method
78         idimension = getdim_relmat(.false.)
79
80         ! Zero & unpack CMO and ZY vectors
81         CALL UPKCMO(CMO,uCMO)
82         IF (NOSIM.GT.0) THEN
83            CALL RSPZYM(NOSIM,BOVECS,BOV)
84            CALL DSCAL(NOSIM*N2ORBX,-1.0d0,BOV,1)
85         END IF
86
87         allocate(mqvec(nosim*idimension))
88         mqvec = 0.0d0
89         allocate(fvvec(nosim*idimension))
90         fvvec = 0.0d0
91
92         ! Determine electric field/potential vector for perturbed
93         ! density matrices
94         call get_fvvec(idimension=idimension,    &
95                        nsim=nosim,   &
96                        udv=udv,      &
97                        cmo=ucmo,     &
98                        work=work,    &
99                        lwork=lwork,  &
100                        fvvec1=fvvec, &
101                        udvtr=udvtr,  &
102                        bovecs=bov)
103
104         ! Allocate and compute Relay matrix
105         idimension = getdim_relmat(.true.)
106         allocate(relmat(idimension))
107         CALL READ_RELMAT(RELMAT)
108         DO I=1,NOSIM
109            IOFF = (I-1)*idimension
110            CALL DGEMV('N',idimension,idimension,1.0d0,RELMAT,idimension,              &
111                       FVVEC(IOFF + 1),1,0.0d0,MQVEC(IOFF + 1),1)
112           if (iprtlvl > 14) then
113              write(lupri, '(/,2x,a,i0)') &
114                  '*** Computed MQ vector start 1st-order density ', i
115              do j = 1, idimension
116                 write(lupri, '(i8, f18.8)') j, MQVEC(IOFF + j)
117              end do
118              write(lupri, '(/,2x,a)') '*** Computed MQ vector end ***'
119           end if
120         END DO
121         deallocate(relmat)
122
123         ! compute xy contributions from induced dipoles and charges
124         if (trplet) then
125            call get_xyvec(ucmo,  &
126                           idimension,  &
127                           nosim, &
128                           rxyt,  &
129                           work,  &
130                           lwork, &
131                           mqvec)
132         else
133            call get_xyvec(ucmo,  &
134                           idimension,  &
135                           nosim, &
136                           rxy,   &
137                           work,  &
138                           lwork, &
139                           mqvec)
140         end if
141      END IF
142
143      ! add qm/np/mm contributions to transformed resp. vectors
144      if (trplet) then
145         ! fixme: i wonder whether this is right, rxy is then never computed
146         call slvsor(.true.,.false.,nosim,udvtr,evecs(1,1),rxy)
147         call slvsor(.true.,.true.,nosim,udv,evecs(1,1),rxyt)
148      else
149         call slvsor(.true.,.true.,nosim,udv,evecs(1,1),rxy)
150      endif
151
152      if (allocated(mqvec)) deallocate(mqvec)
153      if (allocated(fvvec)) deallocate(fvvec)
154      if (allocated(ucmo))  deallocate(ucmo)
155      if (allocated(bov))   deallocate(bov)
156      if (allocated(rxy))   deallocate(rxy)
157      if (allocated(rxyt))  deallocate(rxyt)
158
159   end subroutine
160
161   subroutine qmcmm_qr(vec1,vec2,etrs,xindx,zym1,zym2,              &
162                           udv,work,lwork,kzyvr,kzyv1,kzyv2,            &
163                           igrsym,isymv1,isymv2,cmo,mjwop,              &
164                           ispin0,ispin1,ispin2)
165
166      use qmcmm, only: getdim_relmat, read_relmat
167
168#include "maxorb.h"
169#include "inforb.h"
170#include "infdim.h"
171#include "infinp.h"
172#include "infvar.h"
173#include "infrsp.h"
174#include "infpri.h"
175#include "rspprp.h"
176#include "infcr.h"
177#include "inftap.h"
178#include "qrinf.h"
179#include "mxcent.h"
180#include "priunit.h"
181#include "wrkrsp.h"
182#include "orgcom.h"
183#include "ccinftap.h"
184#include "nuclei.h"
185#include "infpar.h"
186#include "qmnpmm.h"
187
188      integer :: kzyvr
189      integer :: kzyv1
190      integer :: kzyv2
191      integer :: lwork,igrsym,isymv1,isymv2,ispin0,ispin1,ispin2
192      real(8) :: etrs(kzyvr),xindx(*)
193      real(8) :: udv(nashdi,nashdi)
194      real(8) :: zym1(*),zym2(*),work(lwork),cmo(*)
195      real(8) :: vec1(kzyv1),vec2(kzyv2)
196      integer :: mjwop(2,maxwop,8)
197      integer :: i, j, ioff
198      logical   lcon, lorb, lref
199      integer :: idimension
200      integer :: isymt, isymv, isymst, jspin, nsim
201      integer :: nzyvec, nzcvec
202      integer :: ISYMDN, idum
203      real(8) :: ovlap
204
205      real(8), allocatable :: cref(:)
206      real(8), allocatable :: tres(:)
207      real(8), allocatable :: ucmo(:)
208      real(8), allocatable :: tlma(:)
209      real(8), allocatable :: tlmb(:)
210      real(8), allocatable :: trmo(:)
211      real(8), allocatable :: utr(:)
212      real(8), allocatable :: mqvec1(:)
213      real(8), allocatable :: mqvec2(:)
214      real(8), allocatable :: fvvec1(:)
215      real(8), allocatable :: fvvec2(:)
216      real(8), allocatable :: relmat(:)
217
218!     Allocate arrays for response
219      allocate(cref(ncref))
220      allocate(tres(n2orbx))
221      allocate(ucmo(norbt*nbast))
222      allocate(tlma(n2orbx))
223      allocate(tlmb(n2orbx))
224      allocate(trmo(nnorbx))
225      allocate(utr(n2orbx))
226!     Initialize allocated arrays
227      cref = 0.0d0
228      tres = 0.0d0
229      ucmo = 0.0d0
230      tlma = 0.0d0
231      tlmb = 0.0d0
232      trmo = 0.0d0
233      utr = 0.0d0
234!     Reset symmetry variables
235      NSIM  = 1
236      ISYMT = 1
237!     Get the reference state
238      CALL GETREF(CREF,MZCONF(1))
239!     Unpack the response vectors
240      CALL GTZYMT(NSIM,VEC1,KZYV1,ISYMV1,ZYM1,MJWOP)
241      CALL GTZYMT(NSIM,VEC2,KZYV2,ISYMV2,ZYM2,MJWOP)
242!     Unpack symmetry blocked CMO
243      CALL UPKCMO(CMO,UCMO)
244!     Non-iterative method
245      IF (.NOT.MQITER) THEN
246        idimension = getdim_relmat(.false.)
247!       Allocate FV and MQ vectors for all perturbed densities
248        allocate(fvvec1(nsim*idimension))
249        allocate(fvvec2(nsim*idimension))
250
251!       compute fv vectors for second order pertubed densities
252        call get_fvvec(idimension=idimension,     &
253                       nsim=nsim,     &
254                       udv=udv,       &
255                       cmo=ucmo,      &
256                       work=work,     &
257                       lwork=lwork,   &
258                       fvvec1=fvvec1, &
259                       fvvec2=fvvec2, &
260                       isymt=isymt,   &
261                       isymv1=isymv1, &
262                       isymv2=isymv2, &
263                       zym1=zym1,     &
264                       zym2=zym2)
265
266!       Allocate and compute Relay matrix
267        idimension = getdim_relmat(.true.)
268        allocate(relmat(idimension))
269        CALL READ_RELMAT(RELMAT)
270!       Determine induced induced dipoles and charges
271        allocate(mqvec1(nsim*idimension))
272        allocate(mqvec2(nsim*idimension))
273        mqvec1 = 0.0d0
274        mqvec2 = 0.0d0
275        DO I=1,NSIM
276           IOFF = (I-1)*idimension
277           CALL DGEMV('N',idimension,idimension,1.0d0,RELMAT,idimension,              &
278     &                FVVEC1(IOFF + 1),1,0.0d0,MQVEC1(IOFF + 1),1)
279           CALL DGEMV('N',idimension,idimension,1.0d0,RELMAT,idimension,              &
280     &                FVVEC2(IOFF + 1),1,0.0d0,MQVEC2(IOFF + 1),1)
281          if (iprtlvl > 14) then
282             write(lupri, '(/,2x,a,i0)') &
283                 '*** Computed MQ vector start v1 1st-order density ', i
284             do j = 1, idimension
285                write(lupri, '(i8, f18.8)') j, MQVEC1(IOFF+j)
286             end do
287             write(lupri, '(/,2x,a)') '*** Computed MQ vector end ***'
288             write(lupri, '(/,2x,a,i0)') &
289                 '*** Computed MQ vector start v2 1st-order density ', i
290             do j = 1, idimension
291                write(lupri, '(i8, f18.8)') j, MQVEC2(IOFF+j)
292             end do
293             write(lupri, '(/,2x,a)') '*** Computed MQ vector end ***'
294          end if
295        END DO
296        deallocate(mqvec1)
297        deallocate(mqvec2)
298        deallocate(relmat)
299
300        ! determine mm region contribution to qm region potential from
301        ! second order density
302        call get_xyvec(ucmo,   &
303                       idimension,   &
304                       nsim,   &
305                       tres,   &
306                       work,   &
307                       lwork,  &
308                       fvvec1, &
309                       fvvec2, &
310                       isymt,  &
311                       isymv2, &
312                       zym2)
313
314      ELSE
315!       FIX ME: ITERATIVE METHOD
316      END IF
317!     Set up paramterers for quadratic response gradient formation
318      ISYMDN = 1
319      OVLAP  = 1.0d0
320      JSPIN  = 0
321      ISYMV  = IREFSY
322      ISYMST = MULD2H(IGRSYM,IREFSY)
323      IF ( ISYMST .EQ. IREFSY ) THEN
324         LCON = ( MZCONF(IGRSYM) .GT. 1 )
325      ELSE
326         LCON = ( MZCONF(IGRSYM) .GT. 0 )
327      END IF
328      LORB   = ( MZWOPT(IGRSYM) .GT. 0 )
329      LREF = .TRUE.
330      NZYVEC = NCREF
331      NZCVEC = NCREF
332!     Compute gradient
333      CALL RSP1GR(NSIM,KZYVR,IDUM,JSPIN,IGRSYM,JSPIN,ISYMV,ETRS,        &
334     &            CREF,NZYVEC,NZCVEC,OVLAP,ISYMDN,UDV,           &
335     &            TRES,XINDX,MJWOP,WORK,lwork,            &
336     &            LORB,LCON,LREF)
337!
338
339      deallocate(cref)
340      deallocate(tres)
341      deallocate(ucmo)
342      deallocate(tlma)
343      deallocate(tlmb)
344      deallocate(trmo)
345      deallocate(utr)
346
347      end subroutine
348
349
350
351
352
353
354
355   ! computes electric field/potential vector generated by first or second
356   ! order perturbed density matrix
357   subroutine get_fvvec(idimension,   &
358                        nsim,   &
359                        udv,    &
360                        cmo,    &
361                        work,   &
362                        lwork,  &
363                        fvvec1, &
364                        fvvec2, &
365                        udvtr,  &
366                        bovecs, &
367                        isymt,  &
368                        isymv1, &
369                        isymv2, &
370                        zym1,   &
371                        zym2)
372
373      ! size of electric field/potential vector
374      integer, intent(in)              :: idimension
375      ! number of perturbed density matrices
376      integer, intent(in)              :: nsim
377      real(8), intent(in)              :: udv(nashdi, nashdi)
378      real(8), intent(in)              :: cmo(*)
379      real(8), intent(inout)           :: work(*)
380      integer, intent(in)              :: lwork
381      ! electric field/potential vector at np/mm centers
382      real(8), intent(inout)           :: fvvec1(*)
383      real(8), intent(inout), optional :: fvvec2(*)
384      real(8), intent(in),    optional :: udvtr(n2ashx)
385      real(8), intent(in),    optional :: bovecs(*)
386      integer, intent(in),    optional :: isymt
387      integer, intent(in),    optional :: isymv1
388      integer, intent(in),    optional :: isymv2
389      real(8), intent(in),    optional :: zym1(*)
390      real(8), intent(in),    optional :: zym2(*)
391
392#include "dummy.h"
393#include "priunit.h"
394#include "qmnpmm.h"
395#include "inforb.h"
396#include "infdim.h"
397#include "mxcent.h"
398#include "iratdef.h"
399#include "nuclei.h"
400#include "orgcom.h"
401#include "qm3.h"
402#include "infrsp.h"
403
404      logical :: tofile
405      logical :: trimat
406      logical :: exp1vl
407      real(8) :: dipole_origin_save(3)
408      integer :: i
409      integer :: j
410      integer :: ioff
411      integer :: joff
412      integer :: istart
413      integer :: iblk
414      integer :: nocomp
415      integer :: kpatom
416      integer :: isimoff
417      integer :: ixyz
418      real(8) :: f1val, f2val
419      real(8) :: fact
420      real(8) :: tac
421
422      real(8), external :: slvtlm
423      real(8), external :: slvqlm
424
425      real(8), allocatable :: intao(:)
426      real(8), allocatable :: trmo(:)
427      real(8), allocatable :: utr(:)
428      real(8), allocatable :: tlma(:)
429      real(8), allocatable :: tlmb(:)
430      real(8), allocatable :: utrx(:)
431      real(8), allocatable :: urxac(:)
432      integer, allocatable :: intrep(:)
433      integer, allocatable :: intadr(:)
434      character(8), allocatable :: labint(:)
435
436      ! save origin coordinates
437      dipole_origin_save = diporg
438
439      call dzero(fvvec1, idimension*nsim)
440      if (present(fvvec2)) then
441         call dzero(fvvec2, idimension*nsim)
442      end if
443
444      allocate(intao(3*nnbasx))
445      allocate(trmo(nnorbx))
446      allocate(utr(n2orbx))
447      if (present(fvvec2)) then
448         allocate(tlma(n2orbx))
449         allocate(tlmb(n2orbx))
450      else
451         allocate(utrx(n2orbx))
452         allocate(urxac(n2ashx))
453      end if
454      allocate(intrep(9*mxcent))
455      allocate(intadr(9*mxcent))
456      allocate(labint(9*mxcent))
457
458      kpatom = 0
459      tofile = .false.
460      trimat = .true.
461      exp1vl = .false.
462      runqm3 = .true.
463
464      ! compute electric field if needed
465      if (donppol) then
466
467         nocomp = 3
468
469         ! loop over perturbed densities
470         do i = 1, nsim
471            ioff = (i-1)*idimension
472            isimoff = (i-1)*n2orbx+1
473
474            ! loop over np centers
475            do j = 1, tnpatm
476               joff = ioff+(j-1)*3
477               diporg = npcord(:, j)
478
479               intao = 0.0d0
480               call get1in(intao,'NEFIELD',nocomp,work,     &
481                           lwork,labint,intrep,intadr,j,tofile,kpatom,    &
482                           trimat,dummy,exp1vl,dummy,0)
483
484               do ixyz = 1, 3
485
486                  trmo = 0.0d0
487                  utr = 0.0d0
488                  if (.not. present(fvvec2)) then
489                     utrx = 0.0d0
490                     urxac = 0.0d0
491                  end if
492
493                  ! transform integrals
494                  call uthu(intao((ixyz-1)*nnbasx + 1),trmo,cmo,work,nbast,norbt)
495                  call dsptsi(norbt,trmo,utr)
496
497                  if (present(fvvec2)) then
498                     ! determine electric field component size
499                     f1val = 0.0d0
500                     f2val = 0.0d0
501                     if (isymt.eq.isymv1) then
502                        tlma = 0.0d0
503                        call oith1(isymv1,zym1,utr,tlma,isymt)
504                        call melone(tlma,1,udv,1.0d0,f1val,200,'qmnpqro')
505                        fvvec1(joff+ixyz) = f1val
506                     end if
507                     if (isymt.eq.muld2h(isymv1,isymv2)) then
508                        tlma = 0.0d0
509                        tlmb = 0.0d0
510                        call oith1(isymv1,zym1,utr,tlma,isymt)
511                        call oith1(isymv2,zym2,tlma,tlmb,isymv2)
512                        call melone(tlmb,1,udv,1.0d0,f2val,200,'qmnpqro')
513                        fvvec2(joff+ixyz) = f2val
514                     endif
515                  else
516                     call onexh1(bovecs(isimoff),utr,utrx)
517                     if (nasht.gt.0) call getacq(utrx,urxac)
518                     if (trplet) then
519                         fvvec1(joff+ixyz) = slvtlm(udvtr,urxac,utrx,tac)
520                     else
521                         fvvec1(joff+ixyz) = slvqlm(udv,urxac,utrx,tac)
522                     endif
523                  end if
524               end do
525            end do
526         end do
527      end if
528
529      if (donpcap) then
530
531         istart = 0
532         if (donppol) istart = 3*tnpatm
533
534         nocomp = 1
535
536         ! loop over perturbed first order densities
537         do i = 1, nsim
538            ioff = (i-1)*idimension
539            isimoff = (i-1)*n2orbx+1
540
541            ! loop over np centers
542            do j = 1, tnpatm
543               joff = ioff+istart+j
544               diporg = npcord(:, j)
545
546               intao = 0.0d0
547               call get1in(intao,'NPETES ',nocomp,work,     &
548                           lwork,labint,intrep,intadr,j,tofile,kpatom,    &
549                           trimat,dummy,exp1vl,dummy,0)
550
551               trmo = 0.0d0
552               utr = 0.0d0
553               if (.not. present(fvvec2)) then
554                  utrx = 0.0d0
555                  urxac = 0.0d0
556               end if
557
558               ! transform integrals
559               call uthu(intao,trmo,cmo,work,nbast, norbt)
560               call dsptsi(norbt,trmo,utr)
561
562               if (present(fvvec2)) then
563                  f1val = 0.0d0
564                  f2val = 0.0d0
565                  if (isymt.eq.isymv1) then
566                     tlma = 0.0d0
567                     call oith1(isymv1,zym1,utr,tlma,isymt)
568                     call melone(tlma,1,udv,1.0d0,f1val,200,'qmnpqro')
569                     fvvec1(joff) = f1val
570                  end if
571                  if (isymt.eq.muld2h(isymv1,isymv2)) then
572                     tlma = 0.0d0
573                     tlmb = 0.0d0
574                     call oith1(isymv1,zym1,utr,tlma,isymt)
575                     call oith1(isymv2,zym2,tlma,tlmb,isymv2)
576                     call melone(tlmb,1,udv,1.0d0,f2val,200,'qmnpqro')
577                     fvvec2(joff) = f2val
578                  endif
579               else
580                  call onexh1(bovecs(isimoff),utr,utrx)
581                  if (nasht.gt.0) call getacq(utrx,urxac)
582                  if (trplet) then
583                      fvvec1(joff) = slvtlm(udvtr,urxac,utrx,tac)
584                  else
585                      fvvec1(joff) = slvqlm(udv,urxac,utrx,tac)
586                  endif
587               end if
588            end do
589
590            ! set lagrangian for charge equilibration
591            do iblk = 1, tnpblk
592               fvvec1(ioff + idimension) = fvvec1(ioff + idimension) + npchrg(iblk)
593            end do
594            if (present(fvvec2)) then
595               do iblk = 1, tnpblk
596                  fvvec2(ioff + idimension) = fvvec2(ioff + idimension) + npchrg(iblk)
597               end do
598            end if
599         end do
600      end if
601
602      ! print final fv vector
603      do i = 1, nsim
604         if (iprtlvl > 14 .and. .not. mqiter) then
605            ioff = (i - 1)*idimension + 1
606            write(lupri, '(/,2x,a,i0)') &
607                '*** Computed FV vector start 1 ', i
608            do j = 1, idimension
609               write(lupri, '(i8, f18.8)') j, fvvec1(IOFF + j - 1)
610            end do
611            write(lupri, '(/,2x,a)') '*** Computed FV vector end ***'
612            if (present(fvvec2)) then
613               write(lupri, '(/,2x,a,i0)') &
614                   '*** Computed FV vector start 2 ', i
615               do j = 1, idimension
616                  write(lupri, '(i8, f18.8)') j, fvvec2(IOFF + j - 1)
617               end do
618               write(lupri, '(/,2x,a)') '*** Computed FV vector end ***'
619            end if
620         end if
621      end do
622
623      if (allocated(intao))  deallocate(intao)
624      if (allocated(trmo))   deallocate(trmo)
625      if (allocated(utr))    deallocate(utr)
626      if (allocated(tlma))   deallocate(tlma)
627      if (allocated(tlmb))   deallocate(tlmb)
628      if (allocated(utrx))   deallocate(utrx)
629      if (allocated(urxac))  deallocate(urxac)
630      if (allocated(intrep)) deallocate(intrep)
631      if (allocated(intadr)) deallocate(intadr)
632      if (allocated(labint)) deallocate(labint)
633
634      runqm3 = .false.
635
636      ! restore origin coordinates
637      diporg = dipole_origin_save
638
639   end subroutine
640
641
642   ! computes contribution to xy vector from induced dipoles moments and charge
643   subroutine get_xyvec(cmo,     &
644                        idimension,    &
645                        nsim,    &
646                        fvec,    &
647                        work,    &
648                        lwork,   &
649                        fmqvec1, &
650                        fmqvec2, &
651                        isymt,   &
652                        isymv2,  &
653                        zym2)
654
655      ! molecular orbital coefficients
656      real(8), intent(in)           :: cmo(*)
657      ! size of electric field/potential vector
658      integer, intent(in)           :: idimension
659      ! number of perturbed density matrice
660      integer, intent(in)           :: nsim
661      integer, intent(in)           :: lwork
662      real(8), intent(inout)        :: fvec(*)
663      real(8), intent(inout)        :: work(*)
664      ! induced dipole moments and charges
665      real(8), intent(in)           :: fmqvec1(*)
666      real(8), intent(in), optional :: fmqvec2(*)
667      integer, intent(in), optional :: isymt
668      integer, intent(in), optional :: isymv2
669      real(8), intent(in), optional :: zym2(*)
670
671#include "dummy.h"
672#include "qmnpmm.h"
673#include "inforb.h"
674#include "infdim.h"
675#include "mxcent.h"
676#include "iratdef.h"
677#include "nuclei.h"
678#include "orgcom.h"
679#include "qm3.h"
680#include "infrsp.h"
681
682      logical :: tofile
683      logical :: trimat
684      logical :: exp1vl
685      real(8) :: dipole_origin_save(3)
686      integer :: i
687      integer :: j
688      integer :: ioff
689      integer :: joff
690      integer :: istart
691      integer :: kpatom
692      integer :: nocomp
693      integer :: isimoff
694      real(8) :: fact
695      real(8) :: fact1
696      real(8) :: fact2
697      integer :: ixyz
698
699      real(8), allocatable :: utr(:)
700      real(8), allocatable :: trmo(:)
701      real(8), allocatable :: intao(:)
702      real(8), allocatable :: tlma(:)
703      integer, allocatable :: intrep(:)
704      integer, allocatable :: intadr(:)
705      character(8), allocatable :: labint(:)
706
707      ! save origin coordinates
708      dipole_origin_save = diporg
709
710      allocate(utr(n2orbx))
711      allocate(trmo(nnorbx))
712      allocate(intao(3*nnbasx))
713      if (present(fmqvec2)) then
714          allocate(tlma(n2orbx))
715      end if
716      allocate(intrep(9*mxcent))
717      allocate(intadr(9*mxcent))
718      allocate(labint(9*mxcent))
719
720      kpatom = 0
721      tofile = .false.
722      trimat = .true.
723      exp1vl = .false.
724      runqm3 = .true.
725
726      ! induced dipole moment in np region interaction with qm region
727      if (donppol .and. novdamp) then
728         nocomp = 3
729         do i = 1, nsim
730            ioff = (i - 1)*idimension
731            isimoff = (i - 1)*n2orbx + 1
732            do j = 1, tnpatm
733               joff = ioff + (j - 1)*3
734               diporg = npcord(:, j)
735               intao = 0.0d0
736               call get1in(intao, 'NEFIELD', nocomp, work,                   &
737                           lwork, labint, intrep, intadr, j, tofile, kpatom, &
738                           trimat, dummy, exp1vl, dummy, 0)
739
740               ! loop over x, y, z
741               do ixyz = 1, 3
742
743                  ! transform integrals
744                  utr = 0.0d0
745                  trmo = 0.0d0
746                  call uthu(intao((ixyz - 1)*nnbasx + 1), trmo, cmo, work, nbast, norbt)
747                  call dsptsi(norbt, trmo, utr)
748
749                  ! determine mm region contribution
750                  fact1 = -fmqvec1(joff + ixyz)
751                  if (present(fmqvec2)) then
752                     tlma = 0.0d0
753                     call oith1(isymv2, zym2, utr, tlma, isymt)
754                     call daxpy(n2orbx, fact1, tlma, 1, fvec, 1)
755                     fact2 = -0.5d0*fmqvec2(joff + ixyz)
756                     call daxpy(n2orbx, fact2, tlma, 1, fvec, 1)
757                  else
758                     call daxpy(n2orbx, fact1, utr, 1, fvec(isimoff), 1)
759                  end if
760               end do
761            end do
762         end do
763      end if
764
765      ! induced dipole moment in np region interaction with qm region
766      if (donpcap .and. novdamp) then
767         istart = 0
768         if (donppol) istart = 3*tnpatm
769         nocomp = 1
770         do i = 1, nsim
771            ioff = (i - 1)*idimension
772            isimoff = (i - 1)*n2orbx + 1
773            do j = 1, tnpatm
774               joff = ioff + istart + j
775               diporg = npcord(:, j)
776               intao = 0.0d0
777               call get1in(intao, 'NPETES ', nocomp, work,                   &
778                           lwork, labint, intrep, intadr, j, tofile, kpatom, &
779                           trimat, dummy, exp1vl, dummy, 0)
780
781               ! transform integrals
782               utr = 0.0d0
783               trmo = 0.0d0
784               call uthu(intao, trmo, cmo, work, nbast, norbt)
785               call dsptsi(norbt, trmo, utr)
786
787               ! determine mm region contribution
788               fact1 = fmqvec1(joff)
789               if (present(fmqvec2)) then
790                  tlma = 0.0d0
791                  call oith1(isymv2, zym2, utr, tlma, isymt)
792                  call daxpy(n2orbx, fact1, tlma, 1, fvec, 1)
793                  fact2 = 0.5d0*fmqvec2(joff)
794                  call daxpy(n2orbx, fact2, tlma, 1, fvec, 1)
795               else
796                  call daxpy(n2orbx, fact1, utr, 1, fvec(isimoff), 1)
797               end if
798            end do
799         end do
800      end if
801
802      if (allocated(utr))    deallocate(utr)
803      if (allocated(trmo))   deallocate(trmo)
804      if (allocated(intao))  deallocate(intao)
805      if (allocated(tlma))   deallocate(tlma)
806      if (allocated(intrep)) deallocate(intrep)
807      if (allocated(intadr)) deallocate(intadr)
808      if (allocated(labint)) deallocate(labint)
809
810      ! restore origin coordinates
811      diporg = dipole_origin_save
812
813   end subroutine
814
815end module
816