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 cmodel */ 20 SUBROUTINE CMODEL(CMO,GETCMO) 21C 22C Based on TRGEN in Per-Olof Widmark's SCF program. 23C 24C Purpose: delete 3s component in d; 3s,4p in f; etc. 25C If (not GETCMO) then calculate NORB(8) but do not construct CMO. 26C 27C Written 870710-hjaaj 28C revised 89050 -hjaaj: + GETCMO 29C 890703-hjaaj: MOLFTUC4 30C 31#include "implicit.h" 32C 33 DIMENSION CMO(*) 34 LOGICAL GETCMO 35C 36 PARAMETER ( D0 = 0.0D0, DP5= 0.5D0, D1 = 1.0D0, D2 = 2.0D0) 37 PARAMETER ( D3 = 3.0D0, D4 = 4.0D0, D6 = 6.0D0, D8 = 8.0D0) 38 PARAMETER ( D10=10.0D0, D24=24.0D0, D34=34.0D0, D38 = 38.0D0) 39 PARAMETER ( D40=40.0D0, D60=60.0D0, D74=74.0D0, D1270=1270.0D0) 40C 41C Used from common blocks: 42C 43C INFINP : KDEL, TYPE(*) 44C INFORB : NSYM, NORB(8), NBAS(8) 45C PRIUNIT : IPRSTAT, LUSTAT 46C 47#include "maxorb.h" 48#include "priunit.h" 49#include "infinp.h" 50#include "inforb.h" 51#include "infpri.h" 52C 53 CHARACTER*4 MOLFTUC4, CTMPMO(MXCORB) 54C 55 CALL QENTER('CMODEL') 56 IF (IPRSTAT .GE. 99) THEN 57 WRITE (LUSTAT,'(//A/)') ' --- Test output from CMODEL' 58 WRITE (LUSTAT,*) 'KDEL,GETCMO ',KDEL,GETCMO 59 END IF 60C 61 IF (KDEL.EQ.0) THEN 62 IND = 1 63 DO 100 ISYM = 1,NSYM 64 NORB(ISYM) = NBAS(ISYM) 65 IF (GETCMO) CALL DUNIT(CMO(IND),NBAS(ISYM)) 66 IND = IND + NORB(ISYM)*NBAS(ISYM) 67 100 CONTINUE 68 ELSE 69 IF (IPRSIR .GT. 0 .AND. GETCMO) THEN 70 WRITE(LUPRI,*) 71 WRITE(LUPRI,*) 's-component of d-orbitals deleted' 72 WRITE(LUPRI,*) 'p-component of f-orbitals deleted' 73 WRITE(LUPRI,*) 's and d-components of g-orbitals deleted' 74 END IF 75 DO 150 JBAS = 1,NBAST 76 CTMPMO(JBAS) = TYPE(JBAS) 77 150 CONTINUE 78 ITR0 = 0 79 JBAS0 = 0 80 DO 200 ISYM = 1,NSYM 81 NORB(ISYM) = 0 82 DO 210 JBAS = 1,NBAS(ISYM) 83C 84C ... CONSTRUCT PROPER D ORBITALS 85C 86 IF(MOLFTUC4(CTMPMO(JBAS+JBAS0)).EQ.'XX ') THEN 87 NORB(ISYM)=NORB(ISYM)+2 88 JBAS1=JBAS+1 89211 IF(MOLFTUC4(CTMPMO(JBAS1+JBAS0)).NE.'YY ') THEN 90 JBAS1=JBAS1+1 91 IF(JBAS1.GT.NBAS(ISYM)) GOTO 901 92 GOTO 211 93 END IF 94 CTMPMO(JBAS1+JBAS0)='KURT' 95 JBAS2=JBAS1+1 96212 IF(MOLFTUC4(CTMPMO(JBAS2+JBAS0)).NE.'ZZ ') THEN 97 JBAS2=JBAS2+1 98 IF(JBAS2.GT.NBAS(ISYM)) GOTO 901 99 GOTO 212 100 END IF 101 CTMPMO(JBAS2+JBAS0)='KURT' 102 IF(GETCMO) THEN 103 DO 219 I=1,2*NBAS(ISYM) 104 CMO(I+ITR0) = D0 105219 CONTINUE 106 CMO(JBAS +ITR0)= SQRT(DP5) 107 CMO(JBAS1+ITR0)=-SQRT(DP5) 108 ITR0=ITR0+NBAS(ISYM) 109 CMO(JBAS +ITR0)=-D1/SQRT(D6) 110 CMO(JBAS1+ITR0)=-D1/SQRT(D6) 111 CMO(JBAS2+ITR0)= D2/SQRT(D6) 112 ITR0=ITR0+NBAS(ISYM) 113 END IF 114C 115C ... CONSTRUCT PROPER F ORBITALS 116C 117 ELSE IF(MOLFTUC4(CTMPMO(JBAS+JBAS0)).EQ.'XXX ') THEN 118 NORB(ISYM)=NORB(ISYM)+2 119 JBAS1=JBAS+1 120221 IF(MOLFTUC4(CTMPMO(JBAS1+JBAS0)).NE.'XYY ') THEN 121 JBAS1=JBAS1+1 122 IF(JBAS1.GT.NBAS(ISYM)) GOTO 901 123 GOTO 221 124 END IF 125 CTMPMO(JBAS1+JBAS0)='KURT' 126 JBAS2=JBAS1+1 127222 IF(MOLFTUC4(CTMPMO(JBAS2+JBAS0)).NE.'XZZ ') THEN 128 JBAS2=JBAS2+1 129 IF(JBAS2.GT.NBAS(ISYM)) GOTO 901 130 GOTO 222 131 END IF 132 CTMPMO(JBAS2+JBAS0)='KURT' 133 IF(GETCMO) THEN 134 DO 229 I=1,2*NBAS(ISYM) 135 CMO(I+ITR0)=D0 136229 CONTINUE 137 CMO(JBAS +ITR0)= D1/SQRT(D24) 138 CMO(JBAS1+ITR0)=-D3/SQRT(D24) 139 ITR0=ITR0+NBAS(ISYM) 140 CMO(JBAS +ITR0)=-D1/SQRT(D40) 141 CMO(JBAS1+ITR0)=-D1/SQRT(D40) 142 CMO(JBAS2+ITR0)= D4/SQRT(D40) 143 ITR0=ITR0+NBAS(ISYM) 144 END IF 145 ELSE IF(MOLFTUC4(CTMPMO(JBAS+JBAS0)).EQ.'XXY ') THEN 146 NORB(ISYM)=NORB(ISYM)+2 147 JBAS1=JBAS+1 148231 IF(MOLFTUC4(CTMPMO(JBAS1+JBAS0)).NE.'YYY ') THEN 149 JBAS1=JBAS1+1 150 IF(JBAS1.GT.NBAS(ISYM)) GOTO 901 151 GOTO 231 152 END IF 153 CTMPMO(JBAS1+JBAS0)='KURT' 154 JBAS2=JBAS1+1 155232 IF(MOLFTUC4(CTMPMO(JBAS2+JBAS0)).NE.'YZZ ') THEN 156 JBAS2=JBAS2+1 157 IF(JBAS2.GT.NBAS(ISYM)) GOTO 901 158 GOTO 232 159 END IF 160 CTMPMO(JBAS2+JBAS0)='KURT' 161 IF(GETCMO) THEN 162 DO 239 I=1,2*NBAS(ISYM) 163 CMO(I+ITR0)=D0 164239 CONTINUE 165 CMO(JBAS +ITR0)=-D3/SQRT(D24) 166 CMO(JBAS1+ITR0)= D1/SQRT(D24) 167 ITR0=ITR0+NBAS(ISYM) 168 CMO(JBAS +ITR0)=-D1/SQRT(D40) 169 CMO(JBAS1+ITR0)=-D1/SQRT(D40) 170 CMO(JBAS2+ITR0)= D4/SQRT(D40) 171 ITR0=ITR0+NBAS(ISYM) 172 END IF 173 ELSE IF(MOLFTUC4(CTMPMO(JBAS+JBAS0)).EQ.'XXZ ') THEN 174 NORB(ISYM)=NORB(ISYM)+2 175 JBAS1=JBAS+1 176241 IF(MOLFTUC4(CTMPMO(JBAS1+JBAS0)).NE.'YYZ ') THEN 177 JBAS1=JBAS1+1 178 IF(JBAS1.GT.NBAS(ISYM)) GOTO 901 179 GOTO 241 180 END IF 181 CTMPMO(JBAS1+JBAS0)='KURT' 182 JBAS2=JBAS1+1 183242 IF(MOLFTUC4(CTMPMO(JBAS2+JBAS0)).NE.'ZZZ ') THEN 184 JBAS2=JBAS2+1 185 IF(JBAS2.GT.NBAS(ISYM)) GOTO 901 186 GOTO 242 187 END IF 188 CTMPMO(JBAS2+JBAS0)='KURT' 189 IF(GETCMO) THEN 190 DO 249 I=1,2*NBAS(ISYM) 191 CMO(I+ITR0)=D0 192249 CONTINUE 193 CMO(JBAS +ITR0)= DP5 194 CMO(JBAS1+ITR0)=-DP5 195 ITR0=ITR0+NBAS(ISYM) 196 CMO(JBAS +ITR0)=-D3/SQRT(D60) 197 CMO(JBAS1+ITR0)=-D3/SQRT(D60) 198 CMO(JBAS2+ITR0)= D2/SQRT(D60) 199 ITR0=ITR0+NBAS(ISYM) 200 END IF 201C 202C ... CONSTRUCT PROPER G ORBITALS 203C 204 ELSE IF(MOLFTUC4(CTMPMO(JBAS+JBAS0)).EQ.'XXXX') THEN 205 NORB(ISYM)=NORB(ISYM)+3 206 JBAS1=JBAS+1 207251 IF(MOLFTUC4(CTMPMO(JBAS1+JBAS0)).NE.'XXYY') THEN 208 JBAS1=JBAS1+1 209 IF(JBAS1.GT.NBAS(ISYM)) GOTO 901 210 GOTO 251 211 END IF 212 CTMPMO(JBAS1+JBAS0)='KURT' 213 JBAS2=JBAS1+1 214252 IF(MOLFTUC4(CTMPMO(JBAS2+JBAS0)).NE.'XXZZ') THEN 215 JBAS2=JBAS2+1 216 IF(JBAS2.GT.NBAS(ISYM)) GOTO 901 217 GOTO 252 218 END IF 219 CTMPMO(JBAS2+JBAS0)='KURT' 220 JBAS3=JBAS2+1 221253 IF(MOLFTUC4(CTMPMO(JBAS3+JBAS0)).NE.'YYYY') THEN 222 JBAS3=JBAS3+1 223 IF(JBAS3.GT.NBAS(ISYM)) GOTO 901 224 GOTO 253 225 END IF 226 CTMPMO(JBAS3+JBAS0)='KURT' 227 JBAS4=JBAS3+1 228254 IF(MOLFTUC4(CTMPMO(JBAS4+JBAS0)).NE.'YYZZ') THEN 229 JBAS4=JBAS4+1 230 IF(JBAS4.GT.NBAS(ISYM)) GOTO 901 231 GOTO 254 232 END IF 233 CTMPMO(JBAS4+JBAS0)='KURT' 234 JBAS5=JBAS4+1 235255 IF(MOLFTUC4(CTMPMO(JBAS5+JBAS0)).NE.'ZZZZ') THEN 236 JBAS5=JBAS5+1 237 IF(JBAS5.GT.NBAS(ISYM)) GOTO 901 238 GOTO 255 239 END IF 240 CTMPMO(JBAS5+JBAS0)='KURT' 241 IF(GETCMO) THEN 242 DO 259 I=1,3*NBAS(ISYM) 243 CMO(I+ITR0)=D0 244259 CONTINUE 245 CMO(JBAS +ITR0)= -D3/SQRT(D1270) 246 CMO(JBAS1+ITR0)= -D6/SQRT(D1270) 247 CMO(JBAS2+ITR0)= D24/SQRT(D1270) 248 CMO(JBAS3+ITR0)= -D3/SQRT(D1270) 249 CMO(JBAS4+ITR0)= D24/SQRT(D1270) 250 CMO(JBAS5+ITR0)= -D8/SQRT(D1270) 251 ITR0=ITR0+NBAS(ISYM) 252 CMO(JBAS +ITR0)= -D1/SQRT(D74) 253 CMO(JBAS2+ITR0)= D6/SQRT(D74) 254 CMO(JBAS3+ITR0)= D1/SQRT(D74) 255 CMO(JBAS4+ITR0)= -D6/SQRT(D74) 256 ITR0=ITR0+NBAS(ISYM) 257 CMO(JBAS +ITR0)= -D1/SQRT(D38) 258 CMO(JBAS1+ITR0)= D6/SQRT(D38) 259 CMO(JBAS3+ITR0)= -D1/SQRT(D38) 260 ITR0=ITR0+NBAS(ISYM) 261 END IF 262 ELSE IF(MOLFTUC4(CTMPMO(JBAS+JBAS0)).EQ.'XXXY') THEN 263 NORB(ISYM)=NORB(ISYM)+2 264 JBAS1=JBAS+1 265261 IF(MOLFTUC4(CTMPMO(JBAS1+JBAS0)).NE.'XYYY') THEN 266 JBAS1=JBAS1+1 267 IF(JBAS1.GT.NBAS(ISYM)) GOTO 901 268 GOTO 261 269 END IF 270 CTMPMO(JBAS1+JBAS0)='KURT' 271 JBAS2=JBAS1+1 272262 IF(MOLFTUC4(CTMPMO(JBAS2+JBAS0)).NE.'XYZZ') THEN 273 JBAS2=JBAS2+1 274 IF(JBAS2.GT.NBAS(ISYM)) GOTO 901 275 GOTO 262 276 END IF 277 CTMPMO(JBAS2+JBAS0)='KURT' 278 IF(GETCMO) THEN 279 DO 269 I=1,2*NBAS(ISYM) 280 CMO(I+ITR0)=D0 281269 CONTINUE 282 CMO(JBAS +ITR0)= -D1/SQRT(D8) 283 CMO(JBAS1+ITR0)= -D1/SQRT(D8) 284 CMO(JBAS2+ITR0)= D6/SQRT(D8) 285 ITR0=ITR0+NBAS(ISYM) 286 CMO(JBAS +ITR0)= -D1/SQRT(D2) 287 CMO(JBAS1+ITR0)= D1/SQRT(D2) 288 ITR0=ITR0+NBAS(ISYM) 289 END IF 290 ELSE IF(MOLFTUC4(CTMPMO(JBAS+JBAS0)).EQ.'XXXZ') THEN 291 NORB(ISYM)=NORB(ISYM)+2 292 JBAS1=JBAS+1 293271 IF(MOLFTUC4(CTMPMO(JBAS1+JBAS0)).NE.'XYYZ') THEN 294 JBAS1=JBAS1+1 295 IF(JBAS1.GT.NBAS(ISYM)) GOTO 901 296 GOTO 271 297 END IF 298 CTMPMO(JBAS1+JBAS0)='KURT' 299 JBAS2=JBAS1+1 300272 IF(MOLFTUC4(CTMPMO(JBAS2+JBAS0)).NE.'XZZZ') THEN 301 JBAS2=JBAS2+1 302 IF(JBAS2.GT.NBAS(ISYM)) GOTO 901 303 GOTO 272 304 END IF 305 CTMPMO(JBAS2+JBAS0)='KURT' 306 IF(GETCMO) THEN 307 DO 279 I=1,2*NBAS(ISYM) 308 CMO(I+ITR0)=D0 309279 CONTINUE 310 CMO(JBAS +ITR0)= -D3/SQRT(D34) 311 CMO(JBAS1+ITR0)= -D3/SQRT(D34) 312 CMO(JBAS2+ITR0)= D4/SQRT(D34) 313 ITR0=ITR0+NBAS(ISYM) 314 CMO(JBAS +ITR0)= -D1/SQRT(D10) 315 CMO(JBAS1+ITR0)= D3/SQRT(D10) 316 ITR0=ITR0+NBAS(ISYM) 317 END IF 318 ELSE IF(MOLFTUC4(CTMPMO(JBAS+JBAS0)).EQ.'XXYZ') THEN 319 NORB(ISYM)=NORB(ISYM)+2 320 JBAS1=JBAS+1 321281 IF(MOLFTUC4(CTMPMO(JBAS1+JBAS0)).NE.'YYYZ') THEN 322 JBAS1=JBAS1+1 323 IF(JBAS1.GT.NBAS(ISYM)) GOTO 901 324 GOTO 281 325 END IF 326 CTMPMO(JBAS1+JBAS0)='KURT' 327 JBAS2=JBAS1+1 328282 IF(MOLFTUC4(CTMPMO(JBAS2+JBAS0)).NE.'YZZZ') THEN 329 JBAS2=JBAS2+1 330 IF(JBAS2.GT.NBAS(ISYM)) GOTO 901 331 GOTO 282 332 END IF 333 CTMPMO(JBAS2+JBAS0)='KURT' 334 IF(GETCMO) THEN 335 DO 289 I=1,2*NBAS(ISYM) 336 CMO(I+ITR0)=D0 337289 CONTINUE 338 CMO(JBAS +ITR0)= -D3/SQRT(D34) 339 CMO(JBAS1+ITR0)= -D3/SQRT(D34) 340 CMO(JBAS2+ITR0)= D4/SQRT(D34) 341 ITR0=ITR0+NBAS(ISYM) 342 CMO(JBAS +ITR0)= D3/SQRT(D10) 343 CMO(JBAS1+ITR0)= -D1/SQRT(D10) 344 ITR0=ITR0+NBAS(ISYM) 345 END IF 346C 347C ... REST OF BASIS FUNCTIONS 348C ... IGNORE DELETED TYPES 349C 350 ELSE IF(.NOT.(CTMPMO(JBAS+JBAS0).EQ.'KURT')) THEN 351 NORB(ISYM)=NORB(ISYM)+1 352 IF(GETCMO) THEN 353 DO 295 I=1,NBAS(ISYM) 354 CMO(I+ITR0)=D0 355295 CONTINUE 356 CMO(JBAS+ITR0)=D1 357 ITR0=ITR0+NBAS(ISYM) 358 END IF 359 END IF 360210 CONTINUE 361 JBAS0=JBAS0+NBAS(ISYM) 362200 CONTINUE 363 END IF 364 IF(IPRSTAT.GE.9) THEN 365 WRITE(LUSTAT,'(A,8I4)') ' NBAS:',(NBAS(I),I=1,NSYM) 366 WRITE(LUSTAT,'(A,8I4)') ' NORB:',(NORB(I),I=1,NSYM) 367 IF (GETCMO .AND. IPRSTAT.GE.20) THEN 368 WRITE(LUSTAT,*) 'Final CMO in CMODEL' 369 CALL PRORB(CMO,.FALSE.,LUSTAT) 370 END IF 371 END IF 372 CALL QEXIT('CMODEL') 373 RETURN 374901 CONTINUE 375 WRITE(LUERR,*) 'INFINITE LOOP IN CMODEL DETECTED.' 376 WRITE(LUERR,*) 'CONSTRUCTING FOR ',CTMPMO(JBAS+JBAS0) 377 CALL QTRACE(LUERR) 378 CALL QUIT('ERROR CMODEL, INFINITE LOOP (BASIS TYPE NOT FOUND)') 379C ... end of cmodel. 380 END 381C /* Deck molftuc4 */ 382 CHARACTER*4 FUNCTION MOLFTUC4( TEXT ) 383C ( LEFT-ADJUST IN UPPER CASE ) 384C 385C 5-May-1989 -- hjaaj -- force upper case 386C 3-Jul-1989 -- hjaaj -- + left adjust and remove blanks 387C 12-Mar-1991 -- hjaaj -- removed leading p,d,f,g,... used by HERMIT 388C 18-Mar-1993 -- hjaaj -- translate e.g. 'g211' from HERMIT to 'XXYZ' 389C 390 CHARACTER*4 TEXT, TEXTUC 391 CHARACTER*1 CHRUC 392 INTEGER ILCA, ILCZ, UCSHFT, ITEXT 393 LOGICAL FIRST 394C 395 SAVE FIRST, ILCA, ILCZ, UCSHFT 396 DATA FIRST /.TRUE./ 397C 398 IF (FIRST) THEN 399 ILCA = ICHAR('a') 400 ILCZ = ICHAR('z') 401 UCSHFT = ICHAR('A') - ILCA 402 FIRST = .FALSE. 403 END IF 404 TEXTUC = ' ' 405 J = 0 406 IF (TEXT(1:1) .EQ. 'g') THEN 407C ... handle Cartesian g-orbtals from HERMIT 408C (spherical g-orbitals are named '5g-4',...,'5g+4') 409 READ (TEXT,'(1X,3I1)') L,M,N 410 DO 21 I = 1,L 411 J = J + 1 412 TEXTUC(J:J) = 'X' 413 21 CONTINUE 414 DO 22 I = 1,M 415 J = J + 1 416 TEXTUC(J:J) = 'Y' 417 22 CONTINUE 418 DO 23 I = 1,N 419 J = J + 1 420 TEXTUC(J:J) = 'Z' 421 23 CONTINUE 422 ELSE 423 DO 100 I = 1,4 424 IF (TEXT(I:I) .NE. ' ') THEN 425 ITEXT = ICHAR(TEXT(I:I)) 426 IF (ITEXT .GE. ILCA .AND. ITEXT .LE. ILCZ) THEN 427 CHRUC = CHAR( ITEXT + UCSHFT ) 428 ELSE 429 CHRUC = TEXT(I:I) 430 END IF 431C include only 'X', 'Y', or 'Z' in TEXTUC 432 IF (INDEX('XYZ',CHRUC) .NE. 0) THEN 433 J = J + 1 434 TEXTUC(J:J) = CHRUC 435 END IF 436 END IF 437 100 CONTINUE 438 END IF 439 MOLFTUC4 = TEXTUC 440 RETURN 441 END 442C /* Deck delmo */ 443 SUBROUTINE DELMO(CMO,SCRA,LSCRA,THROVL,CMAXMO,DELEMO) 444C 445C Written 18-Jan-198* by Hans Jorgen Aa. Jensen and Hans Agren 446C Revised 26-Aug-1992 by OV+HJAaJ 447C Revised 2004-Feb.2005, Jan.2007 by HJAaJ 448C 449C Purpose: 450C Obtain initial guess for molecular orbitals by 451C diagonalizing the overlap matrix and discarding 452C numerically linar dependent orbitals 453C (defined by eigenvalue of overlap matrix .lt. THROVL). 454C 455C Input: 456C CMO contains initial molecular orbital coefficients 457C (unit matrix, one where 3s in d etc. are deleted, or ?). 458C THROVL: min overlap eigenvalue 459C CMAXMO: max MO coefficient 460C (Generally CMAXMO \approx sqrt(d1/THROVL) 461C but they can be set independently by user. 462C Therefore we require both criteria fulfilled. /hjaaj jan 07) 463C Output: 464C CMO; molecular orbitals 465C DELEMO set true if orbitals have been deleted 466C 467#include "implicit.h" 468#include "priunit.h" 469 DIMENSION CMO(*),SCRA(*) 470 LOGICAL DELEMO 471 PARAMETER (D1=1.0D0) 472C 473C Used from common blocks: 474C INFORB : NNBAST,... 475C INFDIM : NNBASM, NBASMA 476C INFPRI : LUW4 477C CBIREA : LMULBS 478C R12INT : LAUXBS 479C 480#include "inforb.h" 481#include "infdim.h" 482#include "infpri.h" 483#include "cbirea.h" 484#include "r12int.h" 485C 486 CALL QENTER('DELMO ') 487C 488C Core allocation 489C 490 KOVLP = 1 491 KS1T = KOVLP + NNBAST 492 KSCR1 = KS1T + NNBASM 493 KSCR2 = KSCR1 + NBASMA 494 IF (KSCR2+NBASMA .GT. LSCRA) 495 * CALL ERRWRK('DELMO',KSCR2+NBASMA,LSCRA) 496C 497C Read the overlap matrix in AO-basis. 498C 499 CALL RDONEL('OVERLAP ',.TRUE.,SCRA(KOVLP),NNBAST) 500C 501 DELEMO = .FALSE. 502C 503 ICSYM = 1 504 DO 200 ISYM = 1,NSYM 505 NBASI = NBAS(ISYM) 506 ISSYM = KOVLP + IIBAS(ISYM) 507 IF (LAUXBS) THEN 508 NORBI = NORB2(ISYM) 509 ELSE 510 NORB1(ISYM) = NORB(ISYM) 511 NORBI = NORB(ISYM) 512 ICSYM = ICMO(ISYM) + 1 513 END IF 514 IF (NORBI.EQ.0) GO TO 200 515 IF (LMULBS) THEN 516 IF (LAUXBS) THEN 517 IF (.NOT. (R12ECO .OR. R12CBS)) THEN 518C Overlap matrix elements that belong to the 519C primary basis are zeroed (WK/UniKA/04-11-2002). 520C This routine will then delete the corresponding vectors from CMO. 521 MBSYM = ISSYM 522 MBNUL = MBAS1(ISYM) * (MBAS1(ISYM) + 1) / 2 523 CALL DZERO(SCRA(MBSYM),MBNUL) 524 MBSYM = MBSYM + MBNUL 525 DO 202 K = 1, MBAS2(ISYM) 526 CALL DZERO(SCRA(MBSYM),MBAS1(ISYM)) 527 MBSYM = MBSYM + MBAS1(ISYM) + K 528 202 CONTINUE 529 END IF 530 ELSE 531C Overlap matrix elements that belong to the 532C secondary basis are zeroed (WK/UniKA/04-11-2002). 533C This routine will then delete the corresponding vectors from CMO. 534 MBSYM = ISSYM + MBAS1(ISYM) * (MBAS1(ISYM) + 1) / 2 535 MBNUL = NNBAS(ISYM) - MBSYM + ISSYM 536 CALL DZERO(SCRA(MBSYM),MBNUL) 537 ENDIF 538 END IF 539C 540C Update SCR(KS1T) to the new CMO vectors: 541C 542 CALL UTHU(SCRA(ISSYM),SCRA(KS1T),CMO(ICSYM),SCRA(KSCR1), 543 * NBASI,NORBI) 544C 545C Diagonalize overlap matrix ( S * U = U * Seig ) 546C 547 CALL JACO_THR(SCRA(KS1T),CMO(ICSYM),NORBI,NORBI,NBASI,0.0D0) 548C CALL JACO_THR(F,V,NB,NMAX,NROWV,THR_JACO) 549 II = 0 550 DO 175 I = 1,NORBI 551 II = II + I 552 SCRA(KS1T-1+I)=SCRA(KS1T-1+II) 553 175 CONTINUE 554 CALL ORDER2(CMO(ICSYM),SCRA(KS1T),NORBI,NBASI) 555 IF (IPRI6 .GT. 10 .OR. P6FLAG(38)) THEN 556 IF (LMULBS .AND. LAUXBS) THEN 557 WRITE (LUPRI,'(/A,I5/)') 558 & ' DELMO: eigenvalues of overlap matrix for '// 559 & 'auxiliary basis for symmetry',ISYM 560 ELSE 561 WRITE (LUPRI,'(/A,I5/)') 562 & ' DELMO: eigenvalues of overlap matrix for symmetry',ISYM 563 END IF 564 WRITE (LUPRI,'(1P,5D15.5)') (SCRA(KS1T-1+I),I=1,NORBI) 565 END IF 566C 567C Delete orbitals with small ( THROVL ) eigenvalues 568C 569 IDEL = 0 570 DO 275 J = 1, NORBI 571 EIGVAL = SCRA(KS1T + J - 1) 572 JCMO = ICSYM + (J-1)*NBASI 573 ICMOMAX = IDAMAX(NBASI,CMO(JCMO),1) 574 CMOMAX = ABS( CMO(JCMO-1+ICMOMAX) ) / SQRT(EIGVAL) 575 IF (EIGVAL .LT. THROVL .OR. CMOMAX .GT. CMAXMO) THEN 576 IDEL = (NORBI + 1 - J) 577C 578 IF (LAUXBS) THEN 579C ... secondary basis 580 JDEL = IDEL 581 IF (.NOT.DELEMO) WRITE (LUW4,'(//A/A/A,1P,D10.2)') 582 * '@ Auxiliary orbitals are deleted during canonical'// 583 * ' orthonormalization', 584 * '@ because of detected numerical linear dependence.', 585 * '@ Eigenvalue threshold for num. lin. dep. =',THROVL 586 WRITE (LUW4,'(/A,I4,A,I5/A/,(1P,5D10.2))') 587 * '@',IDEL,' MO components are deleted in symmetry',ISYM, 588 * ' Overlap eigenvalues of the deleted orbitals:', 589 * (SCRA(KS1T + J - 2 + I),I = 1,IDEL) 590 DELEMO = .TRUE. 591 IF (ISYM .LT. NSYM) THEN 592 ICTO = ICSYM + (NORB2(ISYM)-IDEL)*NBASI 593 ICFROM= ICSYM + NORB2(ISYM)*NBASI 594 DO 260 JSYM = ISYM+1, NSYM 595 DO 250 I = 0, NORB2(JSYM)*NBAS(JSYM) - 1 596 CMO(ICTO + I) = CMO(ICFROM + I) 597 250 CONTINUE 598 ICTO = ICTO + NORB2(JSYM)*NBAS(JSYM) 599 ICFROM = ICFROM + NORB2(JSYM)*NBAS(JSYM) 600 260 CONTINUE 601 END IF 602 NORB2(ISYM) = NORB2(ISYM) - IDEL 603 ELSE 604C ... primary basis 605 JDEL = IDEL - MBAS2(ISYM) 606 IF (JDEL .GT. 0) THEN 607 IF (.NOT.DELEMO) WRITE (LUW4,'(//A/A/A,1P,D10.2)') 608 * '@ Primary orbitals are deleted during canonical'// 609 * ' orthonormalization', 610 * '@ because of detected numerical linear dependence'// 611 & ' or too large coefficients.', 612 * '@ Eigenvalue threshold for num. lin. dep. =',THROVL 613 WRITE (LUW4,'(/A,I4,A,I5/A/,(1P,5D10.2))') 614 * '@',JDEL,' MO components are deleted in symmetry',ISYM, 615 * ' Overlap eigenvalues of the deleted orbitals:', 616 * (SCRA(KS1T + J - 2 + I),I = 1,JDEL) 617 WRITE (LUW4,'(/A,1P,2D10.2)') 618 & ' Max MO coeff. of last deleted MO & limit:', 619 & CMOMAX,CMAXMO 620 621 DELEMO = .TRUE. 622 END IF 623 IF (ISYM .LT. NSYM) THEN 624 ICTO = ICMO(ISYM) + (NORB(ISYM)-IDEL)*NBASI 625 ICFROM= ICMO(ISYM+1) 626 NCMOVE= NCMOT - ICFROM 627 DO 251 I = 1,NCMOVE 628 CMO(ICTO + I) = CMO(ICFROM + I) 629 251 CONTINUE 630 DO 261 JSYM=ISYM+1,NSYM 631 ICMO(JSYM) = ICMO(JSYM) - IDEL*NBASI 632 261 CONTINUE 633 END IF 634 NORB(ISYM) = NORB(ISYM) - IDEL 635 NCMOT = NCMOT - IDEL*NBASI 636 END IF 637 GOTO 280 638 END IF 639 275 CONTINUE 640 280 CONTINUE 641C 642C Finish canonical orthonormalization by 643C normalization (C = S**(-0.5)*U = U * Seig**(-0.5) ) 644C 645 IF (LAUXBS) THEN 646 NORBI = NORB2(ISYM) 647 ELSE 648 NORBI = NORB(ISYM) 649 NORB1(ISYM) = NORB(ISYM) 650 END IF 651 DO J = 1, NORBI 652 EIGVAL = SCRA(KS1T + J -1) 653 EIGVAL = D1/SQRT(EIGVAL) 654 CALL DSCAL(NBASI,EIGVAL,CMO(ICSYM+(J-1)*NBASI),1) 655 END DO 656C 657 ICSYM = ICSYM + NORB2(ISYM)*NBASI 658 200 CONTINUE 659C 660C *** Make sure primary orbitals are orthonormal /hjaaj-dec-04 661C (there may be numerical round-off errors ...) 662C (ORTHO cannot be used for secondary orbitals as it uses 663C NORB(i) and not NORB2(i). ) 664C 665 IF (.NOT. LAUXBS) THEN 666 KSMAT = 1 667 KWRK = KSMAT + N2BASX 668 LWRK = LSCRA - KWRK 669 CALL ORTHO(CMO,SCRA(KSMAT),SCRA(KWRK),LWRK) 670 END IF 671C 672C *** OUTPUT SECTION 673C 674 IF (LAUXBS) THEN 675 WRITE(LUPRI,'(/A/A)') 676 * ' Canonical primary and secondary basis sets', 677 * ' ISYM NORB1 MBAS1 NORB2 MBAS2' 678 DO ISYM = 1, NSYM 679 WRITE(LUPRI,'(I6,6X,2I6,6X,2I6)') 680 * ISYM, NORB1(ISYM), MBAS1(ISYM), NORB2(ISYM), MBAS2(ISYM) 681 END DO 682 END IF 683C 684C *** end of subroutine DELMO 685C 686 CALL QEXIT('DELMO ') 687 RETURN 688 END 689C /* Deck reord */ 690 SUBROUTINE REORD(CMO,CSCR,IORD) 691C 692C H.AGREN 19.NOV 84 693C 694C Purpose: Reorder mo:s according to IORD(*) array 695C so new_mo(i) = old_mo(iord(i)). 696C 697C Input: CMO : MO:s in normal order 698C 699C Output: CMO : MO:s in IORD order 700C 701#include "implicit.h" 702 DIMENSION CMO(*),CSCR(*),IORD(*) 703C 704C INFORB : NCMOT, ICMO(8), NBAS(8), NORB(8) 705C 706#include "inforb.h" 707C 708 INEW1 = 1 709 DO 10 ISYM = 1,NSYM 710 ICMO1 = ICMO(ISYM)+1 711 NBASI = NBAS(ISYM) 712 IORBI = IORB(ISYM) 713 DO 20 I = 1,NORB(ISYM) 714 IOLD1 = ICMO1 + NBASI*(IORD(IORBI+I) - (IORBI+1)) 715 CALL DCOPY(NBASI,CMO(IOLD1),1,CSCR(INEW1),1) 716 INEW1 = INEW1 + NBASI 717 20 CONTINUE 718 10 CONTINUE 719C 720 CALL DCOPY(NCMOT,CSCR,1,CMO,1) 721C 722 RETURN 723 END 724C /* Deck ortho */ 725 SUBROUTINE ORTHO(CMO,S,SIN,LSIN) 726C 727C Original: CASSCF RELEASE 79 11 23 728C Revisions: 729C 4-Apr-1984 hjaaj 730C Apr-1985 hjaaj (symort, and flag(32) for control) 731C 5-Jul-1989 hjaaj (use PRORB to print orbitals) 732C 733C OBJECTIVE : 734C ORTHOGONALIZES TRIAL MOLECULAR ORBITALS 735C TRANSFERRED FROM SIRINP VIA THE VECTOR CMO 736C 737C SUBROUTINES CALLED: 738C NORM (GRAM-SCHMIDT ORTHONORMALIZATION) 739C SYMORT (SYMMETRICAL ORTHONORMALIZATION) 740C MOLLAB (OVERLAP MATRIX ON LUONEL) 741C 742#include "implicit.h" 743 DIMENSION CMO(*),S(*),SIN(LSIN) 744C 745 PARAMETER (D0 = 0.0D0) 746C 747C Used from common blocks: 748C INFINP : NWARN,CMAXMO,THROVL,? 749C INFORB : NNBAST,NCMOT 750C INFPRI : P4FLAG(*),P6FLAG(*) 751C 752#include "maxorb.h" 753#include "priunit.h" 754#include "infinp.h" 755#include "inforb.h" 756#include "infpri.h" 757C 758 LOGICAL PROVLP 759 SAVE PROVLP 760 DATA PROVLP /.TRUE./ 761C 762 CALL QENTER('ORTHO ') 763C 764C ***** READ OVERLAP MATRIX S FROM LUONEL ***** 765C 766 CALL RDONEL('OVERLAP ',.TRUE.,S,NNBAST) 767C 768 IF (PROVLP .AND. P6FLAG(38)) THEN 769 PROVLP = .FALSE. 770 WRITE(LUPRI,'(/A)') 771 * ' (ORTHO) Overlap between AO basis functions :' 772 CALL OUTPKB(S,NBAS,NSYM,-1,LUPRI) 773 END IF 774C 775C ***** ORTHOGONALIZE ORBITALS SYMMETRY BY SYMMETRY ***** 776C 777 ISTBAS = 0 778 ISTS = 1 779 ISTC = 1 780 DO 100 ISYM=1,NSYM 781 NORBI=NORB(ISYM) 782 NBASI=NBAS(ISYM) 783 IF(NORBI.EQ.0) GO TO 101 784C 785C 786 IF (.NOT.FLAG(32)) THEN 787 CALL NORM(S(ISTS),CMO(ISTC),NBASI,NORBI,SIN,THROVL,IRETUR) 788C ... error exit with IRETUR.ne.0 if norm**2 < THROVL 789C after Gram-Schmidt orthogonalization to prev. vectors 790C hjaaj may 2000: use THROVL instead of fixed THNORM in CALL NORM 791 ELSE 792 CALL SYMORT(CMO(ISTC),S(ISTS),NBASI,NORBI,SIN,LSIN,IRETUR) 793 IF (IRETUR .NE. 0 .AND. IRETUR .NE. 8888) THEN 794C ... if (not converged and not numerical round-off 795C errors) then ... 796 NWARN = NWARN + 1 797 WRITE (LUPRI,4020) ISYM,IRETUR 798 END IF 799 CALL NORM(S(ISTS),CMO(ISTC),NBASI,NORBI,SIN,THROVL,IRETUR) 800C ... 951201: now always call NORM to ensure orthonormality; 801C Cases have been seen where SYMORT have deviation from 802C orthonormality of the order 1.0D-8 because of numerical 803C round-off errors (IRETUR = 8888) 804 END IF 805 IF (IRETUR.NE.0) THEN 806 WRITE (LUPRI,1000) 807 WRITE(LUW4,4010) ISYM, IRETUR 808 IF (LUPRI.NE.LUW4) WRITE(LUPRI,4010) ISYM, IRETUR 809 WRITE(LUERR,4010) ISYM, IRETUR 810 CALL PRORB(CMO,.FALSE.,LUPRI) 811 GO TO 5000 812 END IF 813C 814 101 ISTBAS = ISTBAS + NBASI 815 ISTS = ISTS + NBASI*(NBASI+1)/2 816 ISTC = ISTC + NORBI*NBASI 817 100 CONTINUE 818C 819 4010 FORMAT(/'@*** ORTHO-FATAL ERROR *** Linear dependency in', 820 & ' symmetry =',I3,', CODE =',I4) 821 4020 FORMAT(/'@*** ORTHO-WARNING *** Symmetric orthonormalization'// 822 & ' failed for symmetry',I2, 823 & /'@ Will attempt Gram-Schmidt orthonormalization.') 824 4030 FORMAT(/'@(ORTHO) This error message will now be suppressed', 825 & ' because it has been given max. no. of times (= 3 times).') 826C 827 IMAX = IDAMAX(NCMOT,CMO,1) 828 CMAX = ABS( CMO(IMAX) ) 829c IF (CMAX .GE. CMAXMO) P6FLAG(6) = .TRUE. 830C 831C ***** PRINT MOLECULAR ORBITALS ON UNIT LUPRI ***** 832C 833 IF (P6FLAG(6)) THEN 834 IF (.NOT.FLAG(32)) THEN 835 WRITE (LUPRI,1000) 836 ELSE 837 WRITE (LUPRI,1010) 838 END IF 839 CALL PRORB(CMO,.FALSE.,LUPRI) 840 END IF 841 1000 FORMAT(/' Trial molecular orbitals Gram-Schmidt orthogonalized.') 842 1010 FORMAT(/' Trial molecular orbitals symmetrically orthogonalized.') 843C 844C ***** PRINT MOLECULAR ORBITALS ON UNIT LUW4 ***** 845C 846 IF (P4FLAG(12) .AND. ( LUW4.NE.LUPRI .OR. .NOT.P6FLAG(6) )) THEN 847 IF (.NOT.FLAG(32)) THEN 848 WRITE(LUW4,1000) 849 ELSE 850 WRITE(LUW4,1010) 851 END IF 852 CALL PRORB(CMO,.TRUE.,LUW4) 853 END IF 854C 855 IF (CMAX.GE.CMAXMO) GO TO 104 856C 857 CALL QEXIT('ORTHO ') 858 RETURN 859C 860C ***** ERROR BRANCHES 861C ***** (LINEAR DEPENDENCIES IN ATOMIC BASIS SET) 862C 863C 864 104 CONTINUE 865 DO ISYM = 1,NSYM 866 NORBI = NORB(ISYM) 867 NBASI = NBAS(ISYM) 868 NCMOI = NORBI*NBASI 869 IMAX = IDAMAX(NCMOI,CMO(ICMO(ISYM)+1),1) 870 CMAX = ABS( CMO(ICMO(ISYM)+IMAX) ) 871 IF (CMAX .GT. CMAXMO) THEN 872 IORBI = (IMAX-1)/NBASI + 1 873 WRITE(LUERR,3010) CMAX,ISYM,IORBI,NORBI,CMAXMO 874 WRITE(LUW4 ,3010) CMAX,ISYM,IORBI,NORBI,CMAXMO 875 IF (LUPRI.NE.LUW4) 876 & WRITE(LUPRI,3010) CMAX,ISYM,IORBI,NORBI,CMAXMO 877 END IF 878 END DO 879 WRITE (LUW4,3020) 880 IF (LUPRI.NE.LUW4) WRITE(LUPRI,3020) 881 3010 FORMAT(/' *** ORTHO-FATAL ERROR ***', 882 & /' Largest molecular orbital coefficient/sym',1P,D20.6,3I5, 883 & /' This number is larger than accepted limit',D20.6) 884 3020 FORMAT( 885 & /' Significant loss of accuracy is probable', 886 & /' in transformation of two-electron integrals,' 887 & //'*** Program stops here. *** Your options:', 888 & /' 1) modify basis set,', 889 & /' 2) decrease numerical linear dependency by increasing', 890 & ' ".AO DELETE" limit, or' 891 & /' 3) take the chance and increase', 892 & ' limit with ".CMOMAX" in "*ORBITAL INPUT".') 893 GO TO 5000 894C 895C 896 5000 CALL QTRACE(LUERR) 897 CALL QUIT('*** ERROR *** FATAL ERROR IN ORTHO') 898C 899C end of ORTHO 900C 901 END 902C /* Deck h1mo */ 903 SUBROUTINE H1MO(CMO,SH1,SCRA,LSCRA) 904C 905C Written 15-Apr-1984 by Hans Jorgen Aa. Jensen and Hans Agren 906C Last revision 8-Oct-1984 hjaaj 907C 908C Purpose: 909C Obtain initial guess for molecular orbitals by 910C diagonalizing the one-electron Hamiltonian H1. 911C 912C Input: 913C H1; the one-electron Hamiltonian in AO-basis. 914C 915C Output: 916C CMO; molecular orbitals 917C 918#include "implicit.h" 919 DIMENSION CMO(*),SH1(*),SCRA(*) 920C 921C Used from common blocks: 922C INFINP : FLAG(32) 923C 924#include "maxash.h" 925#include "maxorb.h" 926#include "infinp.h" 927#include "inforb.h" 928#include "infind.h" 929#include "infpri.h" 930#include "infdim.h" 931C 932 LOGICAL LSAVE1,LSAVE2,LSAVE3 933C 934 CALL QENTER('H1MO ') 935C 936C Step 1: Schmidt orthogonalize atomic orbitals 937C (code: flag(32) false) 938C 939C CMO contains initial matrix 940C - either unit matrix or one where some orbitals are deleted 941C (because of e.g. num. lin. dep. or removal of 3s in d). 942C 943 LSAVE1 = P4FLAG(12) 944 LSAVE2 = P6FLAG(6) 945 LSAVE3 = FLAG(32) 946 P4FLAG(12)= .FALSE. 947 P6FLAG(6) = .FALSE. 948 FLAG(32) = .FALSE. 949 CALL ORTHO(CMO,SH1,SCRA,LSCRA) 950C 951C Step 2: Diagonalize one-electron Hamiltonian 952C 953C 954C Get one-electron Hamiltonian 955C 956 CALL SIRH1(SH1,SCRA,LSCRA) 957C 958 KH1T = 1 959 KSCR1 = KH1T + IROW(NBASMA+1) 960 KSCR2 = KSCR1 + NBASMA 961 DO 200 ISYM = 1,NSYM 962 NORBI = NORB(ISYM) 963 IF (NORBI.EQ.0) GO TO 200 964 NBASI = NBAS(ISYM) 965 ISSYM = IIBAS(ISYM) + 1 966 ICSYM = ICMO(ISYM) + 1 967C 968 CALL UTHU(SH1(ISSYM),SCRA(KH1T),CMO(ICSYM),SCRA(KSCR1), 969 * NBASI,NORBI) 970C CALL UTHU(H,HT,U,WRK,NAO,NMO) 971C 972 CALL JACO_THR(SCRA(KH1T),CMO(ICSYM),NORBI,NORBI,NBASI,0.0D0) 973C CALL JACO_THR(F,V,NB,NMAX,NROWV,THR_JACO) 974 II = 0 975 DO 175 I = 1,NORBI 976 II = II + I 977 SCRA(KH1T-1+I)=SCRA(KH1T-1+II) 978 175 CONTINUE 979 CALL ORDER (CMO(ICSYM),SCRA(KH1T),NORBI,NBASI) 980C 981 200 CONTINUE 982C 983C Step 3: Reorthogonalize new mo's using Gram-Schmidt 984C 985 CALL ORTHO(CMO,SH1,SCRA(1),LSCRA) 986 P4FLAG(12)= LSAVE1 987 P6FLAG(6) = LSAVE2 988 FLAG(32) = LSAVE3 989C 990C *** end of subroutine H1MO 991C 992 CALL QEXIT('H1MO ') 993 RETURN 994 END 995C /* Deck h1occ */ 996 SUBROUTINE H1OCC(CMO,WRK,KFRSAV,LFRSAV) 997C 998C Automatic determination of initial HF-occupation 999C from diagonal elements of one-electron Hamiltonian H1. 1000C 1001C Written 25-Aug-1995 by Hans Jorgen Aa. Jensen 1002C 1003C Based in part of modifications originally made in H1MO 1004C by K.Ruud-May 1995 (H1MO now restored to previous content). 1005C 1006C Input: 1007C CMO; molecular orbitals 1008C 1009#include "implicit.h" 1010 DIMENSION CMO(*),WRK(*) 1011C 1012C Used from common blocks: 1013C SCBRHF : IOPRHF 1014C INFORB : NRHF(), NNBAST,NNORBT,NBAST,... 1015C INFIND : IROW() 1016C 1017#include "maxash.h" 1018#include "maxorb.h" 1019#include "priunit.h" 1020#include "scbrhf.h" 1021#include "inforb.h" 1022#include "infind.h" 1023#include "infpri.h" 1024C 1025 CALL QENTER('H1OCC ') 1026C 1027 KFREE = KFRSAV 1028 LFREE = LFRSAV 1029 CALL MEMGET('REAL',KHMO,NNORBT,WRK,KFREE,LFREE) 1030 CALL MEMGET('REAL',KHAO,NNBAST,WRK,KFREE,LFREE) 1031C 1032C ***** retrieve atomic ONE ELECTRON HAMILTONIAN matrix 1033C 1034 CALL SIRH1(WRK(KHAO),WRK(KFREE),LFREE) 1035C 1036C Transform ONE ELECTRON HAMILTONIAN to MO basis 1037C 1038 CALL MEMGET('REAL',KWRK,N2BASX,WRK,KFREE,LFREE) 1039 CALL UTHUB(WRK(KHAO),WRK(KHMO),CMO,WRK(KWRK),NSYM,NBAS,NORB) 1040C 1041C ***** Test output of ONE ELECTRON HAMILTONIAN matrices ***** 1042C 1043 IF (IPRI6 .GE. 15) THEN 1044 WRITE(LUPRI,1000) 1045 CALL OUTPKB(WRK(KHAO),NBAS,NSYM,-1,LUPRI) 1046 WRITE(LUPRI,1200) 1047 CALL OUTPKB(WRK(KHMO),NORB,NSYM,-1,LUPRI) 1048 END IF 1049 CALL MEMREL('H1OCC.1',WRK,KFRSAV,KHAO,KFREE,LFREE) 1050C 1051 1000 FORMAT(/' H1OCC: TEST OF ONE ELECTRON HAMILTONIAN (AO-basis)') 1052 1200 FORMAT(/' H1OCC: TEST OF ONE ELECTRON HAMILTONIAN (MO-basis)') 1053C 1054C Extract H1 diagonal in WRK(KH1D) and 1055C orbital symmetries in WRK(KSMO) 1056C 1057 CALL MEMGET('REAL',KH1D,NORBT,WRK,KFREE,LFREE) 1058 CALL MEMGET('REAL',KSMO,NORBT,WRK,KFREE,LFREE) 1059 DO 200 ISYM = 1,NSYM 1060 NORBI = NORB(ISYM) 1061 IF (NORBI.EQ.0) GO TO 200 1062 JHMO = KHMO-1 + IIORB(ISYM) 1063 JH1D = KH1D-1 + IORB(ISYM) 1064 JSMO = KSMO-1 + IORB(ISYM) 1065 DO 175 I = 1,NORBI 1066 WRK(JH1D+I) = WRK(JHMO+IROW(I+1)) 1067 WRK(JSMO+I) = ISYM 1068 175 CONTINUE 1069 200 CONTINUE 1070C 1071C Sort according to diagonal element and 1072C determine HF occupation. 1073C 1074 CALL ORDER(WRK(KSMO),WRK(KH1D),NORBT,1) 1075 CALL IZERO(NRHF,8) 1076 MOCC = NRHFEL/2 1077 DO 20 I = 1, MOCC 1078 ISYM = NINT(WRK(KSMO-1+I)) 1079 NRHF(ISYM) = NRHF(ISYM) + 1 1080 20 CONTINUE 1081 IF (2*MOCC .NE. NRHFEL) IOPRHF = NINT(WRK(KSMO+MOCC)) 1082C 1083C *** end of subroutine H1OCC 1084C 1085 CALL QEXIT('H1OCC ') 1086 RETURN 1087 END 1088C /* Deck fcmo */ 1089 SUBROUTINE FCMO(MXFCK,CMO,FC,SCRA,LSCRA) 1090C 1091C Written 4-May-1984 by Hans Jorgen Aa. Jensen 1092C Revisions: 1093C 8-Oct-1984 hjaaj 1094C 1-Nov-1984 hjaaj (use NORB(*), not NBAS(*) for FC) 1095C 1985 hjaaj (use NRHF(*) for RHF occupation) 1096C 4-Mar-1997 tsaue include screening 1097C 1098C Purpose: 1099C Do MXFCK restricted Roothaan-Hartree-Fock iterations. 1100C 1101C -- idea: grand-canonical Hartree-Fock can be specified 1102C by NASH(*) = NRHF(*), NASHT = NRHFT, 1103C DV(ij) = delta(ij) occ(i); but NISHT = 0 and NISH(*) = 0. 1104C Then GC Fock matrix is FC + FV = h1 + FV. (860117) 1105C 1106C Input: 1107C CMO; initial molecular orbitals used to build Fock matrix, 1108C assumed to be orthonormal. 1109C MXFCK; maximum number of Fock iterations (always one iteration) 1110C 1111C Output: 1112C CMO; molecular orbitals diagonalizing Fock matrix 1113C 1114C Scratch: 1115C FC; the inactive Fock matrix and scratch area for overlap matrix 1116C SCRA; general scratch area 1117C 1118#include "implicit.h" 1119 DIMENSION CMO(*),FC(*),SCRA(*) 1120C 1121C 1122 PARAMETER (D0=0.0D0, EMYCNV = 1.D-4) 1123#include "dummy.h" 1124C 1125C Used from common blocks: 1126C INFINP : 1127C INFOPT : EPOT,EMY,EMCSCF 1128C SCBRHF : NFRRHF(*) 1129C INFIND : ...,ISSMO(*),? 1130C 1131#include "maxash.h" 1132#include "maxorb.h" 1133#include "priunit.h" 1134#include "infinp.h" 1135#include "infopt.h" 1136#include "inforb.h" 1137#include "scbrhf.h" 1138#include "infind.h" 1139#include "infpri.h" 1140#include "gnrinf.h" 1141C 1142 LOGICAL LSAVE4,LSAVE6 1143C 1144 CALL QENTER('FCMO ') 1145 WRITE (LUPRI,'(A//A/A//A,I5/A,8I5)') '1', 1146 * ' Restricted Roothaan-Hartree-Fock iterations', 1147 * ' -------------------------------------------', 1148 * ' Number of electrons :',2*NRHFT, 1149 * ' Orbital occupations :',(NRHF(I),I=1,NSYM) 1150C 1151 NFRHFT = ISUM(NSYM,NFRRHF,1) 1152 IF (NFRHFT .GT. 0) THEN 1153 NWARN = NWARN + 1 1154 WRITE (LUPRI,'(//A/A/)') 1155 * '@ WARNING from FCMO: ".FROZEN" is not implemented'// 1156 * ' for Fock iterations,', 1157 * '@ Fock iterations abandoned.' 1158 GO TO 9999 1159 END IF 1160 IF (NASHT .GT. 0) THEN 1161 NWARN = NWARN + 1 1162 WRITE (LUPRI,'(//A/A/)') 1163 * '@ WARNING from FCMO: HSROHF and ROHF are not implemented'// 1164 * ' for Fock iterations,', 1165 * '@ Fock iterations abandoned.' 1166 GO TO 9999 1167 END IF 1168C 1169 LSAVE4 = P4FLAG(12) 1170 LSAVE6 = P6FLAG(6) 1171 P4FLAG(12)= .FALSE. 1172 P6FLAG(6) = .FALSE. 1173 DO 5 ISYM = 1,NSYM 1174 ISWAP = NRHF(ISYM) 1175 NRHF(ISYM) = NISH(ISYM) 1176 NISH(ISYM) = ISWAP 1177 5 CONTINUE 1178 ISWAP = NRHFT 1179 NRHFT = NISHT 1180 NISHT = ISWAP 1181 IEXIT = 0 1182 ITFCK = 0 1183 EMY = D0 1184 100 CONTINUE 1185 ITFCK = ITFCK + 1 1186 EMYOLD = EMY 1187C 1188C Step 1: Construct inactive Fock matrix 1189C 1190 CALL FCKMAT(.TRUE.,DUMMY,CMO,EMY,FC,DUMMY,SCRA,LSCRA) 1191C CALL FCKMAT(ONLYFC,DV,CMO,EMY,FC,FV,WRK,LFREE) 1192 IF (SUPSYM) CALL AVEFCK(FC) 1193 EMCSCF = EPOT + EMY 1194 IF (IEXIT .EQ. 1) GO TO 500 1195 IF (IPRSIR .GT. 0) WRITE (LUPRI,1730) ITFCK,EMY,EMCSCF 1196C 1197C Step 2: Diagonalize inactive Fock-matrix: 1198C 1199 DO 200 ISYM = 1,NSYM 1200 NORBI = NORB(ISYM) 1201 IF (NORBI.EQ.0) GO TO 200 1202 IORBI = IORB(ISYM) 1203 NBASI = NBAS(ISYM) 1204 ISSYM = IIORB(ISYM) + 1 1205 ICSYM = ICMO(ISYM) + 1 1206C 1207 CALL JACO_THR(FC(ISSYM),CMO(ICSYM),NORBI,NORBI,NBASI,0.0D0) 1208C CALL JACO_THR(F,VEC,NB,NMAX,NROWV,THR_JACO) 1209C 1210 II = ISSYM - 1 1211 DO 175 I=1,NORBI 1212 II = II + I 1213 SCRA(I)=FC(II) 1214 175 CONTINUE 1215 CALL ORDRSS(CMO(ICSYM),SCRA,ISSMO(IORBI+1),NORBI,NBASI) 1216 IF (IPRSIR .GE. 5) THEN 1217 IF (IPRSIR .GE. 12) THEN 1218 IEND = NORBI 1219 ELSE 1220 IEND = MIN(NORBI,NISH(ISYM)+2) 1221 END IF 1222 WRITE (LUPRI,1740) ISYM 1223 WRITE (LUPRI,1750) (SCRA(I),I=1,IEND) 1224 END IF 1225 200 CONTINUE 1226 IF (SUPSYM) CALL AVEORD() 1227C ... remake ISSORD() as ISSMO() may have changed in ORDRSS 1228C 1229 1730 FORMAT(//' Fock iteration number',I3, 1230 * '; inactive energy:',F25.12, 1231 * /T28,'total energy:',F25.12) 1232 1740 FORMAT(/' Fock eigenvalues symmetry',I2/) 1233 1750 FORMAT(4X,5F15.8) 1234C 1235C Step 3: Reorthogonalize new mo's 1236C 1237 CALL ORTHO(CMO,FC,SCRA(1),LSCRA) 1238C 1239C Another Fock iteration? 1240C 1241ckr IF (ITFCK .GT. 1 .AND. EMY .GT. EMYOLD) THEN 1242C 1243C Stop because of oscillation, i.e. energy has gone up. 1244C One more iteration has been taken. This should lead to 1245C lower energy if the oscillation is typical. 1246C 1247ckr WRITE (LUPRI,'(//2A/)') ' *** Energy is oscillating,', 1248ckr * ' exit from Roothaan-Hartree-Fock iterations.' 1249ckr ELSE 1250 IF (ITFCK .GT. 1 .AND. abs(EMYOLD-EMY) .LE. EMYCNV) THEN 1251 WRITE (LUPRI,'(//2A,1P,D10.2/)') 1252 * ' *** Roothaan-Hartree-Fock energy difference', 1253 * ' converged to',(EMYOLD-EMY) 1254 ELSE IF (ITFCK .LT. MXFCK) THEN 1255 GO TO 100 1256C ^----------------- 1257 ELSE 1258 WRITE (LUPRI,'(//A/)') 1259 & ' *** Max. number of iterations reached.' 1260 END IF 1261ckr END IF 1262 IEXIT = 1 1263 GO TO 100 1264C ... go calculate final energy 1265C 1266C 1267C 1268 500 CONTINUE 1269 WRITE (LUPRI,1730) ITFCK,EMY,EMCSCF 1270 P4FLAG(12)= LSAVE4 1271 P6FLAG(6) = LSAVE6 1272 DO 800 ISYM = 1,NSYM 1273 ISWAP = NRHF(ISYM) 1274 NRHF(ISYM) = NISH(ISYM) 1275 NISH(ISYM) = ISWAP 1276 800 CONTINUE 1277 ISWAP = NRHFT 1278 NRHFT = NISHT 1279 NISHT = ISWAP 1280C 1281C *** end of subroutine FCMO 1282C 1283 9999 CALL QEXIT('FCMO ') 1284 RETURN 1285 END 1286C /* Deck fcvirt */ 1287 SUBROUTINE FCVIRT(CMO,WRK,LFREE) 1288C 1289C 2-Oct-1986 Poul Joergensen 1290C Revised 28-Aug-1995 hjaaj 1291C 1292C Purpose : To use one electron hamiltonian or modified FC 1293C to modify virtual Hartree-Fock orbitals for CI/MC 1294C 1295C Reference for .FC MVO: C.W. Bauschlicher, JCP 72 (1980) 880. 1296C 1297#include "implicit.h" 1298C 1299 DIMENSION CMO(*), WRK(LFREE) 1300C 1301 PARAMETER ( D0 = 0.0D0 ) 1302#include "dummy.h" 1303C 1304C INFORB : NSYM, NNBAST, NNORBT, ... 1305C SCBRHF : IOPRHF, NMVO(), NMVOT 1306C INFIND : IROW(*) 1307C INFPRI : IPRI6 1308C INFIND : SUPSYM 1309C 1310#include "maxash.h" 1311#include "maxorb.h" 1312#include "priunit.h" 1313#include "inforb.h" 1314#include "scbrhf.h" 1315#include "infind.h" 1316#include "infpri.h" 1317#include "infinp.h" 1318C 1319 DIMENSION MISH(8), MRHF(8) 1320C 1321 CALL QENTER('FCVIRT') 1322C 1323 CALL ICOPY(8,NRHF,1,MRHF,1) 1324 IF (IOPRHF .GT. 0) MRHF(IOPRHF) = MRHF(IOPRHF) + 1 1325C ... if NRHF() was not defined in *RHF CALC 1326C it will contain NISH() from *WAVE FUNC 1327C 1328 IF (NMVOT .EQ. 0) THEN 1329 WRITE (LUPRI,'(//A/A)') 1330 & ' --- Modified virtual orbitals generated by diagonalization', 1331 & ' --- of virtual block of one-electron Hamiltonian.' 1332 ELSE 1333 WRITE (LUPRI,'(//A/A/A,8I4)') 1334 & ' --- Modified virtual orbitals generated by diagonalization', 1335 & ' --- of virtual block of core Fock matrix.', 1336 & ' Number of core orbitals in each symmetry : ', 1337 & (NMVO(I),I=1,NSYM) 1338 END IF 1339 WRITE (LUPRI,'(A,8I4/)') 1340 & ' Number of occupied orbitals in each symmetry: ', 1341 & (MRHF(I),I=1,NSYM) 1342C 1343 KFC_MVO = 1 1344 KWRK = KFC_MVO + NNORBT 1345 LNEED= KWRK + 2*NBAST 1346 LWRK = LFREE - KWRK 1347 IF (LNEED .GT. LFREE) CALL ERRWRK('FCVIRT',LNEED,LFREE) 1348C 1349C ***** Calculate a Fock matrix "FC_MVO" based on NMVO(*) occupied orbitals 1350C instead fo NISH(*). 1351C 1352 CALL ICOPY(8,NISH,1,MISH,1) 1353 CALL ICOPY(8,NMVO,1,NISH,1) 1354 MISHT = NISHT 1355 NISHT = NMVOT 1356 CALL FCKMAT(.TRUE.,DUMMY,CMO,EMCMY,WRK(KFC_MVO),DUMMY, 1357 & WRK(KWRK),LWRK) 1358C CALL FCKMAT(ONLYFC,DV,CMO,EMCMY,FC,FV,WRK,LFRSAV) 1359 CALL ICOPY(8,MISH,1,NISH,1) 1360 NISHT = MISHT 1361C 1362C ***** Zero all elements in FC_MVO matrix except virtual Hartree-Fock block 1363C 1364 DO 200 ISYM = 1,NSYM 1365 NORBI = NORB(ISYM) 1366 IF (NORBI.EQ.0) GO TO 200 1367 IORBI = IORB(ISYM) 1368 NBASI = NBAS(ISYM) 1369 NOCCI = MRHF(ISYM) 1370 JXFC = KFC_MVO + IIORB(ISYM) 1371 JCMO = 1 + ICMO(ISYM) 1372C 1373C 1374 DO 167 I = 1 , NORBI 1375 JSTA = JXFC + IROW(I) 1376 JEND = JSTA - 1 + MIN(I,NOCCI) 1377 DO 164 J = JSTA,JEND 1378 WRK(J) = D0 1379 164 CONTINUE 1380 167 CONTINUE 1381C 1382C 1383 CALL JACO_THR(WRK(JXFC),CMO(JCMO),NORBI,NORBI,NBASI,0.0D0) 1384C CALL JACO_THR(F,VEC,NB,NMAX,NROWV,THR_JACO) 1385C 1386 DO 175 I = 1,NORBI 1387 WRK(I) = WRK(JXFC-1 + IROW(I+1)) 1388 175 CONTINUE 1389C 1390 NSSHI = NORBI - NOCCI 1391 IF (NSSHI .GT. 0) THEN 1392 JCMO = JCMO + NOCCI*NBASI 1393C order virtual HARTREE-FOCK orbitals 1394 CALL ORDRSS(CMO(JCMO),WRK(1+NOCCI), 1395 & ISSMO(IORBI+NOCCI+1),NSSHI,NBASI) 1396 END IF 1397C 1398 IF (IPRI6 .GE. 10) 1399 * WRITE (LUPRI,'(/A/A,I3,//,(5(I3,F12.6)))') 1400 * ' Super symmetry and eigenvalues of virtual one-electron'// 1401 * ' Hamiltonian',' Symmetry',ISYM, 1402 * (ISSMO(IORBI+I),WRK(I),I=(NOCCI+1),NORBI) 1403 200 CONTINUE 1404 IF (SUPSYM) CALL AVEORD() 1405C ... remake ISSORD() as ISSMO() may have changed in ORDRSS 1406C 1407 CALL QEXIT('FCVIRT') 1408 RETURN 1409 END 1410C /* Deck prorb */ 1411 SUBROUTINE PRORB(CMO,PROCC,IOUT) 1412C 1413C Written 11-Dec-1983 by ha/hjaaj 1414C Revised 8-Jul-1992 hjaaj 1415C 1416C Purpose: 1417C Print molecular orbital coefficients on unit IOUT. 1418C 1419C Input: 1420C CMO: MO orbital coefficients (symmetry blocked) 1421C PROCC: print only occupied orbitals (unless CMOPRI true) 1422C IOUT: output file unit 1423C 1424#include "implicit.h" 1425 DIMENSION CMO(*) 1426 LOGICAL PROCC 1427C 1428C Used from common blocks: 1429C pgroup : REP 1430C INFINP : CENT,TYPE,SUPSYM,? 1431C INFORB : NSYM,... 1432C INFIND : ISSMO(), 1433C INFPRI : CMOPRI 1434C CBIREA : LMULBS 1435C R12INT : MBAS1(:) 1436C 1437#include "maxorb.h" 1438#include "maxash.h" 1439#include "pgroup.h" 1440#include "infinp.h" 1441#include "inforb.h" 1442#include "infind.h" 1443#include "infpri.h" 1444#include "cbirea.h" 1445#include "r12int.h" 1446C 1447C *** FORMAT statements 1448C 1449 2400 FORMAT(/' Molecular orbitals for symmetry species',I2,' (',A,')' 1450 & /' ------------------------------------------------') 1451 2600 FORMAT(/' Orbital ',7I9) 1452 2620 FORMAT( ' Sup.sym.',7I9) 1453 2800 FORMAT(I4,1X,A4,':',A4,7F9.4) 1454C 1455C *** Print molecular orbitals 1456C 1457 IF (PROCC .AND. .NOT. CMOPRI) THEN 1458 THRPRI = 1.0D-2 1459 WRITE(IOUT,'(/A,F7.4,A)') 1460 & ' (Only coefficients >',THRPRI,' are printed.)' 1461 ELSE 1462 THRPRI = 1.0D-4 1463 END IF 1464C 1465 ISTBAS = 0 1466 DO 400 ISYM = 1,NSYM 1467 NBASI = NBAS(ISYM) 1468 MBASI = NBAS(ISYM) 1469 IF (LMULBS) MBASI = MBAS1(ISYM) 1470C 1471 IF (PROCC .AND. .NOT. CMOPRI) THEN 1472 IF (NBAST .LE. 50) THEN 1473C print all occupied and 10 lowest unoccupied this symmetry 1474 IEND = 0 1475 ELSE 1476C print 10 highest doubly occupied, all active, and 10 lowest unoccupied this symmetry 1477C (output becomes too big otherwise ... /hjaaj Dec08) 1478 IEND = MAX(NISH(ISYM)-10, 0) 1479 END IF 1480 NENDI = MIN(NOCC(ISYM)+10, NORB(ISYM)) 1481 ELSE 1482 IEND = 0 1483 NENDI = NORB(ISYM) 1484 END IF 1485 IF (NENDI.EQ.0) GO TO 300 1486 WRITE(IOUT,2400) ISYM,REP(ISYM-1) 1487C 1488 ICMOI = ICMO(ISYM) 1489 ISTORB = IORB(ISYM) 1490 100 IST =IEND+1 1491 ISTMO =IEND*NBASI+ICMOI 1492 IEND =MIN(IEND+7,NENDI) 1493 IEMO =NBASI*(IEND-1)+ICMOI 1494 WRITE(IOUT,2600) (I,I=IST,IEND) 1495 IF (SUPSYM) WRITE(IOUT,2620) (ISSMO(ISTORB+I),I=IST,IEND) 1496 DO 200 I=1,MBASI 1497 JSMO=ISTMO+I 1498 JEMO=IEMO+I 1499 JJ = 0 1500 DO J=JSMO,JEMO,NBASI 1501 IF ( ABS(CMO(J)) .GE. THRPRI ) JJ = 1 1502 END DO 1503 IF (JJ .EQ. 1) WRITE(IOUT,2800) I,CENT(I+ISTBAS), 1504 * TYPE(I+ISTBAS),(CMO(J),J=JSMO,JEMO,NBASI) 1505 200 CONTINUE 1506 IF (IEND.NE.NENDI) GO TO 100 1507C 1508 300 CONTINUE 1509 ISTBAS = ISTBAS + NBASI 1510 400 CONTINUE 1511C 1512C *** End of subroutine PRORB 1513C 1514 RETURN 1515 END 1516C /* Deck punmo */ 1517 SUBROUTINE PUNMO(IPCTL,CMO,OCC) 1518C 1519C 25-May-1985 hjaaj 1520C l.r. 910715-hjaaj (added formats 7010-7050) 1521C 1522C Purpose: 1523C punch mo coefficients and, conditionally, occupation numbers 1524C to LUPNCH 1525C 1526C Control: 1527C IPCTL .eq. 1: also punch occupation numbers 1528C else : do not punch occupation numbers 1529C 1530#include "implicit.h" 1531#include "dummy.h" 1532 DIMENSION CMO(*),OCC(*) 1533 PARAMETER (D0 = 0.0D0) 1534C 1535C Used from common blocks: 1536C INFORB : NSYM,NORB(*),NBAS(*) 1537C 1538#include "inforb.h" 1539C 1540 CHARACTER*8 LBLDAT(2) 1541C 1542 CALL QENTER('PUNMO') 1543 CALL GETDAT(LBLDAT(1),LBLDAT(2)) 1544 LUPNCH = -1 1545 CALL GPOPEN(LUPNCH,'DALTON.MOPUN','UNKNOWN',' ','FORMATTED', 1546 & IDUMMY,.FALSE.) 1547 REWIND LUPNCH 1548 IF (IPCTL.EQ.1) THEN 1549 WRITE (LUPNCH,'(A,A8,2X,A8,A)') 1550 * '**NATORB ( punched by SIRIUS ',LBLDAT,' )' 1551 ELSE 1552 WRITE (LUPNCH,'(A,A8,2X,A8,A)') 1553 * '**MOLORB ( punched by SIRIUS ',LBLDAT,' )' 1554 END IF 1555C 1556 IEND =0 1557 DO 700 ISYM=1,NSYM 1558 NORBI=NORB(ISYM) 1559 IF (NORBI.EQ.0) GO TO 700 1560 NBASI=NBAS(ISYM) 1561 DO NI=1,NORBI 1562 IST=IEND+1 1563 IEND=IEND+NBASI 1564 IMOMXV = IDAMAX(NBASI,CMO(IST),1) 1565 CMOMXV = ABS(CMO(IST-1+IMOMXV)) 1566 IF (CMOMXV .LT. 1.0D2) THEN 1567 WRITE (LUPNCH,7000) (CMO(NP),NP=IST,IEND) 1568 ELSE IF (CMOMXV .LT. 1.0D3) THEN 1569 WRITE (LUPNCH,7010) (CMO(NP),NP=IST,IEND) 1570 ELSE IF (CMOMXV .LT. 1.0D4) THEN 1571 WRITE (LUPNCH,7020) (CMO(NP),NP=IST,IEND) 1572 ELSE IF (CMOMXV .LT. 1.0D5) THEN 1573 WRITE (LUPNCH,7030) (CMO(NP),NP=IST,IEND) 1574 ELSE IF (CMOMXV .LT. 1.0D6) THEN 1575 WRITE (LUPNCH,7040) (CMO(NP),NP=IST,IEND) 1576 ELSE 1577 WRITE (LUPNCH,7050) (CMO(NP),NP=IST,IEND) 1578 END IF 1579 END DO 1580 NDELI = NBASI-NORBI 1581 IF (NDELI.EQ.0) GO TO 700 1582 DO NI=1,NDELI 1583 WRITE (LUPNCH,7000) (D0,NP=1,NBASI) 1584 END DO 1585 700 CONTINUE 1586C 1587 IF (IPCTL.EQ.1) THEN 1588 WRITE (LUPNCH,'(A,A8,2X,A8,A)') 1589 * '**NATOCC ( punched by SIRIUS ',LBLDAT,' )' 1590 WRITE(LUPNCH,7000) (OCC(NP),NP=1,NORBT) 1591 END IF 1592C 1593 CALL GPCLOSE(LUPNCH,'KEEP') 1594 CALL QEXIT('PUNMO') 1595 RETURN 1596C 1597 7000 FORMAT(4F18.14) 1598 7010 FORMAT(4F18.13) 1599 7020 FORMAT(4F18.12) 1600 7030 FORMAT(4F18.11) 1601 7040 FORMAT(4F18.10) 1602 7050 FORMAT(4G18.11) 1603C 1604C end of PUNMO 1605C 1606 END 1607 1608 SUBROUTINE SIR_VIRTRUNC(WORK,LWORK) 1609! 1610! 7-Nov-2017 Hans Joergen Aa. Jensen 1611! 1612! Remove virtual orbitals according to .VIRTRUNC input 1613! 1614#include "implicit.h" 1615#include "priunit.h" 1616 REAL*8 WORK(LWORK) 1617! 1618! gnrinf.h : WRINDX 1619! infinp.h : THR_VIRTRUNC 1620! inforb.h : NCMOT, ... 1621! inforb.h : NCONF 1622! inftap.h : LUIT1 1623#include "maxorb.h" 1624#include "gnrinf.h" 1625#include "infinp.h" 1626#include "inforb.h" 1627#include "infvar.h" 1628#include "inftap.h" 1629 LOGICAL ANTIS, OLDWOP, FNDLB2 1630 INTEGER NORB_NEW(8), NORBT_NEW, IROW, I 1631 CHARACTER*8 RN_LABEL, STAR8, RTNLBL(2) 1632 1633 IROW(I) = (I*(I-1))/2 1634 1635 KFREE = 1 1636 LFREE = LWORK 1637! retrieve CMO 1638 CALL MEMGET2('REAL','CMO',KCMO,NCMOT,WORK,KFREE,LFREE) 1639 JRDMO = 9 1640 CALL READMO(WORK(KCMO),JRDMO) 1641! calculate R**n integrals in MO basis over secondary/virtual indices 1642! (symmetry packed) 1643 CALL MEMGET2('REAL','RN VIRPK ',KRN_VIR,NNORBT,WORK,KFREE,LFREE) ! NNSSHT has not been calculated, but NNORBT .gt. NNSSHT 1644 CALL MEMGET2('REAL','RN MO PK ',KRN_MO,NNORBT,WORK,KFREE,LFREE) 1645 CALL MEMGET2('REAL','RN AO PK ',KRN_AO,NNBAST,WORK,KFREE,LFREE) 1646 CALL MEMGET2('REAL','RN AO INT',KRNINT,NNBASX,WORK,KFREE,LFREE) 1647 IF (N_in_RN .eq. 2) THEN 1648 RN_LABEL = 'R2 ' 1649 THR_RN = THR_VIRTRUNC ** 2 1650 ELSE IF (N_in_RN .eq. 4) THEN 1651 RN_LABEL = 'R4 ' 1652 THR_RN = THR_VIRTRUNC ** 4 1653 ELSE 1654 WRITE(LUPRI,'(//A,I0)') 1655 & '(SIR_VIRTRUNC) Illegal N_in_RN from .TRUNCVIRT :',N_in_RN 1656 CALL QUIT('SIR_VIRTRUNC: Illegal N_in_RN from .TRUNCVIRT') 1657 END IF 1658 CALL RDPROP(RN_LABEL,WORK(KRNINT),ANTIS) 1659 CALL PKSYM1(WORK(KRNINT),WORK(KRN_AO),NBAS,NSYM,1) 1660 CALL UTHUB(WORK(KRN_AO),WORK(KRN_MO),WORK(KCMO), 1661 & WORK(KFREE),NSYM,NBAS,NORB) 1662 CALL GETVIR(WORK(KRN_MO),WORK(KRN_VIR)) 1663 CALL MEMREL('RN INT',WORK,1,KRN_MO,KFREE,LFREE) 1664 1665! extract virtual (secondary) block in each symmetry, 1666! diagonalize and remove unwanted MO:s 1667 NORB_NEW(:) = NORB(:) 1668 NORBT_NEW = NORBT 1669 IISSH_I = 0 1670 DO ISYM = 1, NSYM 1671 NSSH_I = NSSH(ISYM) 1672 IF (NSSH_I .EQ. 0) CYCLE 1673 NBAS_I = NBAS(ISYM) 1674 JRN_VIR = KRN_VIR + IISSH_I 1675 JCMO = KCMO + ICMO(ISYM) + NBAS_I*NOCC(ISYM) 1676 1677 CALL JACO_THR(WORK(JRN_VIR),WORK(JCMO),NSSH_I,NSSH_I,NBAS_I, 1678 & 1.0D-12) 1679 1680 JCMO_1 = JCMO 1681 JCMO_2 = JCMO 1682 DO I = 1, NSSH_I 1683 II = JRN_VIR-1 + IROW(I+1) 1684 IF (WORK(II) .LT. THR_RN) THEN 1685 IF (JCMO_1 .ne. JCMO_2) THEN 1686 CALL DCOPY(NBAS_I,WORK(JCMO_1),1,WORK(JCMO_2),1) 1687 END IF 1688 JCMO_2 = JCMO_2 + NBAS_I 1689 ELSE ! remove this virtual/secondary orbital 1690 NORB_NEW(ISYM) = NORB_NEW(ISYM) - 1 1691 NORBT_NEW = NORBT_NEW - 1 1692 END IF 1693 JCMO_1 = JCMO_1 + NBAS_I 1694 END DO 1695 1696 IISSH_I = IISSH_I + (NSSH_I*(NSSH_I+1))/2 1697 END DO 1698 1699 IF (NORBT_NEW .EQ. NORBT) THEN 1700 WRITE(LUPRI,'(/A)') 1701 &' .VIRTRUNC: no orbitals removed, SIRIUS.RST not updated.' 1702 GO TO 9000 ! if no orbitals removed, then exit 1703 END IF 1704 1705 WRITE(LUPRI,'(/A,I0,A)') '.VIRTRUNC: ',NORBT-NORBT_NEW, 1706 &' virtual/secondary orbitals removed, SIRIUS.RST updated.' 1707 WRITE(LUPRI,'(A,8I5)') 1708 &' Number of MOs per symmetry reduced from',NORB(1:NSYM) 1709 WRITE(LUPRI,'(A,8I5)') 1710 &' to ',NORB_NEW(1:NSYM) 1711 1712 CALL MEMGET2('REAL','CMO2',KCMO2,NCMOT,WORK,KFREE,LFREE) 1713 JCMO2 = KCMO2 1714 DO ISYM = 1, NSYM 1715 JCMO1 = KCMO + ICMO(ISYM) 1716 NCMO_I = NORB_NEW(ISYM)*NBAS(ISYM) 1717 CALL DCOPY(NCMO_I,WORK(JCMO1),1,WORK(JCMO2),1) 1718 JCMO2 = JCMO2 + NCMO_I 1719 END DO 1720 NORB(:) = NORB_NEW(:) 1721 NORBT = NORBT_NEW 1722 1723 NRS = 0 1724 IF (NCONF .GT. 1) THEN 1725 REWIND LUIT1 1726 IF (FNDLB2('STARTVEC',RTNLBL,LUIT1)) THEN 1727 READ (RTNLBL(1),'(2I4)') NRS, I 1728 CALL MEMGET2('REAL','CVECS',KCVECS,NRS*NCONF, 1729 & WORK,KFREE,LFREE) 1730 DO II = 1, NRS 1731 JCVECS = KCVECS + (II-1)*NCONF 1732 CALL READT(LUIT1,NCONF,WORK(JCVECS)) 1733 END DO 1734 END IF 1735 END IF 1736 1737! update orbital information with the reduced number of molecular orbitals 1738 WRINDX = .TRUE. 1739 OLDWOP = .FALSE. 1740 CALL SIRSET(WORK(KFREE),LFREE,OLDWOP) 1741 IAVERR = 0 1742 CALL AVECHK(IAVERR) 1743 IF (IAVERR .NE. 0) CALL QUIT( 1744 & 'SIR_VIRTRUNC error: inconsistency in sup.sym. averaging') 1745 ! write new basis size info on LUIT1 1746 CALL NEWIT1 1747C 1748C save NRS CI vectors read above (if any): 1749C 1750 IF (NRS .GT. 0) THEN 1751 WRITE (LUIT1) '********',RTNLBL(1),'VIRTRUNCSTARTVEC' 1752 DO II = 1, NRS 1753 JCVECS = KCVECS + (II-1)*NCONF 1754 CALL WRITT(LUIT1,NCONF,WORK(JCVECS)) 1755 END DO 1756 END IF 1757 ! write new reduced set of MO coefficients 1758 WRITE (LUIT1) '******** VIRTRUNCNEWORB ' 1759 CALL WRITT(LUIT1,NCMOT,WORK(KCMO2)) 1760 WRITE (LUIT1) '******** VIRTRUNCEODATA ' 1761 FLUSH(LUIT1) 1762 ! hjaaj 9-Nov-2017: for some strange reason did gfortran 5.5 not empty the output buffer before the QUIT without the flush statement ... 1763 REWIND (LUIT1) 1764 1765 9000 CONTINUE 1766 CALL MEMREL('VIRTRUNC',WORK,1,1,KFREE,LFREE) 1767 1768 CALL QUIT('DALTON: Planned exit because .VIRTRUNC finished.') 1769 1770 RETURN 1771 END 1772C --- end of sirius/sircmo.F --- 1773