1!
2!  Dalton, a molecular electronic structure program
3!  Copyright (C) by the authors of Dalton.
4!
5!  This program is free software; you can redistribute it and/or
6!  modify it under the terms of the GNU Lesser General Public
7!  License version 2.1 as published by the Free Software Foundation.
8!
9!  This program is distributed in the hope that it will be useful,
10!  but WITHOUT ANY WARRANTY; without even the implied warranty of
11!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
12!  Lesser General Public License for more details.
13!
14!  If a copy of the GNU LGPL v2.1 was not distributed with this
15!  code, you can obtain one at https://www.gnu.org/licenses/old-licenses/lgpl-2.1.en.html.
16!
17!
18C
19C  /* Deck cc_molden */
20      subroutine cc_molden_nto(work,lwork)
21C
22C     R. Faber et al., 2018
23C
24C     based on CC_MOLDEN by
25C     Rolf H. Myhre and H. Koch
26C
27C     Purpose: Calculate natural transition orbitals and write them in
28C     molden format to file
29C
30      implicit none
31
32      character(len=*), parameter :: myname = 'CC_MOLDEN_NTO'
33!
34#include "priunit.h"
35#include "ccorb.h"
36#include "ccsdinp.h"
37#include "ccsdsym.h"
38#include "ccfop.h"
39#include "ccsections.h"
40#include "ccexgr.h"
41#include "ccexci.h"
42#include "ccexcinf.h"
43#include "inftap.h"
44
45!
46      integer :: kcmot, knto, kucmo, koccu, krbtr, korbv
47      integer :: kvec1, kumat, kvmat, ksval, kpos
48      integer :: isyma, isymi, isym
49      integer :: excisym, iexci
50      integer :: mindim
51      integer :: nto_count(8), ncount, nfound, ncount2
52      integer :: i_info, idummy, ilist, ilistoff, iopt
53      integer ludena
54      integer :: icmo(8), iorbs(8)
55      integer :: kend1, kend2, kend3
56      integer :: lwrk1, lwrk2, lwrk3
57      integer, intent(in) :: lwork
58!
59      character(len=20):: filename
60      character(len=10):: model
61      character(len=3), parameter :: list_type = 'RE '!'LE '
62!
63#if defined (SYS_CRAY)
64      real, intent(out):: work(lwork)
65      real :: one, half, ddot
66#else
67      double precision, intent(out):: work(lwork)
68      double precision, parameter :: one = 1.d0, half = 0.5d0,
69     &                               zero = 0d0,
70     &                               thresh_nto = 1.0D-2
71      double precision :: ddot, DNRM2
72      double precision trace, dummy
73#endif
74!
75      model ='CCSD      '
76!
77      nto_count = 0
78!
79      ncount  = 0
80      ncount2 = 0
81      do isym = 1, nsym
82         icmo(isym) = ncount
83         iorbs(isym) = ncount2
84         !ncount = ncount + nbas(isym)*norb(isym)
85         ncount = ncount + nbas(isym)*norbs(isym)
86         ncount2 = ncount2 + norbs(isym)
87      end do
88!
89!     Dynamic allocation
90!     ------------------
91!
92      kcmot = 1
93      knto = kcmot + nlamds
94      kucmo = knto + nlamds
95      krbtr = kucmo + nbast*norbts
96      korbv = krbtr + nbast*nbast
97      koccu = korbv + norbts
98      kend1 = koccu + norbts
99      lwrk1 = lwork - kend1 + 1
100!      write(lupri,'(6I5)') kcmot, kucmo, krbtr, korbv, koccu, kend1
101!
102      if (lwrk1 .lt. 0) then
103         stop 'Insufficient work space in '//myname
104      endif
105
106!
107!     Close molden_MOS, we want new files
108!     -----------------------------------
109      call gpclose(LUMOLDEN_MOS,'KEEP')
110!
111!
112!     Read in MO coefficients
113!     -------------------------
114      IF (LUSIFC .LE. 0) CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ',
115     &                               'UNFORMATTED',IDUMMY,.FALSE.)
116      REWIND LUSIFC
117C
118      CALL MOLLAB('TRCCINT ',LUSIFC,LUERR)
119      READ (LUSIFC)
120C
121      READ (LUSIFC)
122      READ (LUSIFC) (WORK(KCMOT+I-1),I=1,nlamds)
123C
124      CALL GPCLOSE(LUSIFC,'KEEP')
125C
126C     Reorder and delete frozen orbitals
127      CALL CMO_REORDER(WORK(KCMOT), WORK(KEND1), LWRK1)
128!
129!     Natural transition orbitals
130!     ---------------------
131!
132      do excisym = 1, nsym
133
134         ! symmetry specific entries
135         kvec1 = kend1
136         kend2 = kvec1 + nt1am(excisym)
137         lwrk2 = lwork - kend2
138
139         ilistoff = ISYOFE(excisym)
140
141         do iexci = 1, nccexci(excisym,1) + nccexci(excisym,3)
142            iopt = 1 ! Read only singles part of vector
143            ilist = ilistoff + iexci
144            ! Read in actual vector
145            call cc_rdrsp(list_type, ilist, excisym, iopt,
146     &                     model, work(kvec1), dummy)
147
148            ! Zero the nto memory
149            call dzero(work(knto), nlamds)
150            call dzero(work(koccu), norbts)
151
152            ! Loop over symmetry blocks of the vector
153            do isymi = 1, nsym
154               isyma = muld2h(isymi,excisym)
155
156               mindim = min(nvir(isyma),nrhf(isymi))
157               ! Skip if there is no work to do
158               if (mindim.lt.1) cycle
159
160               ! Allocate memory
161               kumat = kend2
162               kvmat = kumat + nvir(isyma)**2
163               ksval = kvmat + nrhf(isymi)**2
164               kend3 = ksval + mindim
165               lwrk3 = lwork - kend3
166
167               kpos = kvec1 + it1am(isyma,isymi)
168
169               ! Do the SVD... Maybe we can use 'S' instead of 'A'
170               CALL DGESVD(
171     &            'A','A',nvir(isyma),nrhf(isymi),
172     &            WORK(kpos),
173     &            nvir(isyma), work(ksval),
174     &            work(kumat), nvir(isyma),
175     &            work(kvmat), nrhf(isymi),
176     &            work(kend3), lwrk3, i_info
177     &         )
178               if (i_info .ne. 0) then
179                  write(LUPRI,'(A,i5)') 'ERROR in DGSVD ', i_info
180                  cycle
181               endif
182
183
184               ! work(kumat) now has the left/virtual eigenvectors on
185               ! the columns
186               ! work(kvmat) now has the right/occupied eigevectors on
187               ! the rows
188
189               ! Find the number of NTOs to write
190               nfound = 0
191               do n = 1, mindim
192                  if (work(ksval-1+n) .lt. thresh_nto) exit
193                  nfound = nfound + 1
194               end do
195               ! Ensure that we don't write too many, i.e. more than
196               ! half the number ofbitals
197               nfound = min(nfound, norbs(isymi)/2, norbs(isyma)/2)
198
199               ! Transform the hole NTOs to AO basis
200               call dgemm('N','T', nbas(isymi), nfound, nrhf(isymi),
201     &                    one,
202     &                    work(kcmot+ilmrhf(isymi)), nbas(isymi), ! AO x occ-MO matrix
203     &                    work(kvmat), nrhf(isymi), ! NTO x occ-MO matrix
204     &                    zero, work(knto+icmo(isymi)), nbas(isymi))
205
206               ! Transform the particle NTOs to AO basis
207               kpos = knto + icmo(isyma)
208     &              + nbas(isyma)*(norbs(isyma)-nfound)
209               call dgemm('N','N', nbas(isyma), nfound, nvir(isyma),
210     &                    one,
211     &                    work(kcmot+ilmvir(isyma)), nbas(isyma),
212     &                    work(kumat), nvir(isyma),
213     &                    zero, work(kpos), nbas(isyma))
214
215               ! Put the hole NTO eigenvalues in the occupation vector
216               kpos = koccu + iorbs(isymi)
217               work(kpos:kpos+nfound-1) = - work(ksval:ksval+nfound-1)
218               kpos = koccu + iorbs(isyma) + (norbs(isyma) - nfound)
219               work(kpos:kpos+nfound-1) = work(ksval:ksval+nfound-1)
220
221            end do
222!
223!           Create filename
224!           ---------------
225!
226            write(filename,"(A5,I1,A,I0,A7)")
227     &         "exci_", excisym, '_', iexci, ".molden"
228!
229!           Call molden_mos to write MOs in molden format
230!           ---------------------------------------------
231            call gpopen(LUMOLDEN_MOS,trim(filename),'NEW',' ',
232     &                  'FORMATTED',idummy,.false.)
233            rewind lumolden_mos
234            call molden_mos(1, work(knto), work(koccu), work(krbtr),
235     &                      work(kucmo),work(korbv))
236            call gpclose(LUMOLDEN_MOS,'KEEP')
237         end do
238      end do
239!
240!     open original molden_mos.tmp
241!     ----------------------------
242      call gpopen(LUMOLDEN_MOS,'molden_MOS.tmp','UNKNOWN',' ',
243     &            'FORMATTED',idummy,.false.)
244!
245      return
246      end
247