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