C C $Id: sproc.F,v 1.2 1998/07/16 16:40:46 jjv5 Exp arjan $ C C------------------------------------------------------------------------ C SUBSYSTEM PROCESSING ROUTINES. C C SPROC1 SHOULD BE CALLED WHENEVER A NEW SUBSYSTEM ATOM LIST IS C CREATED. IT SHOULD BE CALLED BEFORE ATTEMPTING TO CREATE OR C UPDATE THE PAIRLIST. SPROC1 REQUIRES THE SUBSYSTEM INFORMATION C NSUB, IATOM1, IATOMS, IABUFF. C C SPROC2 SHOULD BE CALLED AFTER EACH PAIRLIST UPDATE. IT REQUIRES C PAIRLIST ARRAYS IP1 AND IPAIR AND THE SUBSYSTEM ARRAYS ISUB1 C AND ISUBS. C C C****************************************************************************** C C SUBROUTINE SPROC1(IERROR) C C DOES A VARIETY OF SUBSYSTEM SETUP WORK PRIOR TO CREATING C PAIRLIST. C C BUG FIX: 6/29/98 -- S. DIXON. WHEN DUMMY ATOMS APPEAR AFTER THE C FIRST ATOM, THE IS1ATM POINTER WASN'T BEING UPDATED PROPERLY. C THIS OCCURRED IN THE "DO 500 IATM=1,NATOMS LOOP". NOW FIXED. C 7/26/99 FIX ALSO APPLIES TO SPARKLES C IMPLICIT DOUBLE PRECISION (A-H,O-Z) #include "divcon.dim" #include "divcon.h" C C LOCAL: C DIMENSION ISFLAG(MAXSUB) C C C ASSIGN ATOMIC ORBITAL POINTERS FOR EACH SUBSYSTEM: C C IORBPT(K) = GLOBAL POINTER TO FIRST ATOMIC ORBITAL IN SUBSYSTEM K. C IORB1(I) = LOCAL POINTER IN A GIVEN SUBSYSTEM FOR THE FIRST C ATOMIC ORBITAL ON ATOM IATOM(I). C C SET ERROR FLAG IF GLOBAL OR LOCAL STORAGE LIMITS ARE EXCEEDED. C IORBK = 1 NSORB = 0 DO 100 K=1,NSUB IORBPT(K) = IORBK I1 = IATOM1(K) I2 = IATOM1(K+1)-1 NORBS = 0 IORBI = 1 DO 20 I=I1,I2 IATM = IATOMS(I) NORBSI = NATORB(IATNUM(IATM)) NORBS = NORBS + NORBSI IORB1(I) = IORBI IORBI = IORBI + NORBSI 20 CONTINUE NSORB = MAX(NSORB,NORBS) IF(NORBS.GT.MSORB)THEN WRITE(IOUT,'(" SUBSYSTEM ",I4," HAS ",I4," ATOMIC ORBITALS,", . " WHICH EXCEEDS MAXIMUM ALLOWED")') K,NORBS ENDIF IORBK = IORBK + NORBS 100 CONTINUE IORBPT(NSUB+1) = IORBK IF(NSORB.GT.MSORB)THEN IERROR = 1 WRITE(IOUT,'(/" LOCAL STORAGE FOR SUBSYSTEM ATOMIC ORBITALS", . " EXCEEDED --"/" INCREASE MSORB PARAMETER IN", . " divcon.dim TO AT LEAST ",I4)') NSORB ENDIF IF(IORBK.GT.MSVAL)THEN IERROR = 1 WRITE(IOUT,'(/" GLOBAL STORAGE FOR SUBSYSTEM EIGENVALUES", . " EXCEEDED --"/" INCREASE MSVAL PARAMETER IN", . " divcon.dim TO AT LEAST ",I7)') IORBK ENDIF IF(IERROR.NE.0) RETURN C C CREATE LISTS OF SUBSYSTEMS ISUBS IN WHICH EACH ATOM APPEARS AS C A NON-BUFFER ATOM. C IPT = 1 DO 200 IATM=1,NATOMS C C ISUB1(IATM) IS A POINTER TO THE POSITION IN ISUBS THAT C CONTAINS THE FIRST SUBSYSTEM IN WHICH IATM APPEARS AS A C NON-BUFFER ATOM. C ISUB1(IATM) = IPT C C DON'T BOTHER TO LOOK FOR SUBSYSTEMS IF IATM IS NOT A REAL ATOM. C IF(NATORB(IATNUM(IATM)).EQ.0) GO TO 200 DO 180 KSUB=1,NSUB J1 = IATOM1(KSUB) J2 = IATOM1(KSUB+1)-1 C C SEE IF IATM COINCIDES WITH ANY NON-BUFFER ATOMS IN KSUB. C DO 160 J=J1,J2 JATM = IATOMS(J) C C IF JATM IS A HIGHER ATOM NUMBER THAN IATM, WE DON'T NEED TO C CONTINUE THE SEARCH IN KSUB BECAUSE THE ATOM LIST IS SORTED. C IF(JATM.GT.IATM) GO TO 180 C C SKIP CURRENT ATOM IF IT'S A BUFFER IN KSUB. C IF(IABUFF(J).NE.0) GO TO 160 C IF(JATM.EQ.IATM)THEN C C STORE KSUB IN THE LIST FOR IATM AND CHECK THE NEXT SUBSYSTEM. C ISUBS(IPT) = KSUB IPT = IPT + 1 GO TO 180 ENDIF 160 CONTINUE 180 CONTINUE 200 CONTINUE C C STORE OVERFLOW POINTER. C ISUB1(NATOMS+1) = IPT C C NOW CHECK TO MAKE SURE THAT EACH REAL ATOM APPEARS AT LEAST ONCE C AS A NON-BUFFER ATOM. C DO 300 I=1,NATOMS IF(NATORB(IATNUM(I)).EQ.0) GO TO 300 NSUBS = ISUB1(I+1)-ISUB1(I) IF(NSUBS.LE.0)THEN IERROR = 1 WRITE(IOUT,'(" ATOM ",I5," DOES NOT APPEAR AS A NON-BUFFER", . " ATOM IN ANY SUBSYSTEM")') I ENDIF 300 CONTINUE C C FOR EACH ATOM, GENERATE THE LIST OF SUBSYSTEMS IASUBS IN WHICH C THE ATOM APPEARS. IASUBS DIFFERS FROM ISUBS IN THAT IASUBS C DOES NOT DISTINGUISH WHETHER OR NOT THE ATOM HAS BUFFER STATUS. C IPT = 1 DO 500 IATM=1,NATOMS C C POINTER FOR FIRST SUBSYSTEM BELONGING TO IATM: C IS1ATM(IATM) = IPT C C SKIP IATM IF IT'S NOT A REAL ATOM. THIS IS THE BUG FIX. THE GO TO C STATEMENT APPEARED BEFORE THE POINTER ASSIGNMENT IN THE PREVIOUS C VERSION, AND IT CAUSED AN ERROR WITH DUMMY ATOMS -- S. DIXON. C IF(NATORB(IATNUM(IATM)).EQ.0) GO TO 500 C DO 480 KSUB=1,NSUB J1 = IATOM1(KSUB) J2 = IATOM1(KSUB+1)-1 C C SEE IF IATM MATCHES ANY ATOMS IN KSUB. C DO 460 J=J1,J2 JATM = IATOMS(J) C C BECAUSE THE SUBSYSTEM ATOM LISTS ARE SORTED, WE CAN DISCONTINUE C SEARCH IN KSUB IF WE'VE 'PASSED UP' IATM. C IF(JATM.GT.IATM) GO TO 480 C IF(JATM.EQ.IATM)THEN C C IATM IS IN KSUB. STORE KSUB IN ITS LIST. C IASUBS(IPT) = KSUB IPT = IPT + 1 GO TO 480 ENDIF 460 CONTINUE 480 CONTINUE 500 CONTINUE C C OVERFLOW POINTER: C IS1ATM(NATOMS+1) = IPT RETURN END C C C SUBROUTINE SPROC2 C C FOR EACH PAIR OF ATOMS (IATM,JATM) IN THE PAIRLIST, THIS ROUTINE C ASSIGNS THE NUMBER OF SUBSYSTEMS NSHARE IN WHICH BOTH ATOMS APPEAR: C C K = PAIRLIST ENTRY FOR (IATM,JATM). C NSHARE(1,K) = NUMBER OF TIMES (IATM,JATM) APPEAR WITH AT LEAST C ONE ATOM HAVING NON-BUFFER STATUS, AND THE OTHER C ATOM HAVING NON-BUFFER STATUS OR INNER BUFFER LAYER C STATUS. C NSHARE(2,K) = NUMBER OF TIMES (IATM,JATM) APPEAR WITH ANY STATUS. C C C NSHARE(1,K) KEEPS TRACK OF HOW MANY TIMES A PAIR OF ATOMS CONTRIBUTE C TO THE DENSITY MATRIX, WHILE NSHARE(2,K) KEEPS TRACK OF CONTRIBUTIONS C TO THE 1-ELECTRON AND FOCK MATRICES. C C SHOULD BE CALLED AFTER EACH PAIRLIST UPDATE, BUT ONLY AFTER C CALLING IJMAKE. C IMPLICIT DOUBLE PRECISION (A-H,O-Z) #include "divcon.dim" #include "divcon.h" C IF(NATOMS.EQ.1) RETURN C NPAIRS = IP1(NATOMS+1)-1 DO 100 I=1,NPAIRS NSHARE(1,I) = 0 NSHARE(2,I) = 0 100 CONTINUE C C LOOP OVER EACH PAIR OF ATOMS IN EACH SUBSYSTEM AND ACCUMULATE C NSHARE BASED ON BUFFER STATUS. C DO 500 K=1,NSUB NATMS = IATOM1(K+1) - IATOM1(K) IF(NATMS.EQ.1) GO TO 500 IATM0 = IATOM1(K)-1 DO 400 I=2,NATMS I0I = IATM0+I IATM = IATOMS(I0I) IBUFFI = IABUFF(I0I) DO 300 J=1,I-1 I0J = IATM0+J JATM = IATOMS(I0J) C C SKIP THE IATM,JATM PAIR IF IT HASN'T BEEN STORED AS A C BONDED PAIR. C CALL IJFIND(NPAIRS,IATM,JATM,IJADDR) IF(IJADDR.EQ.0) GO TO 300 C IBUFFJ = IABUFF(I0J) IJSUM = IBUFFI + IBUFFJ IF(IJSUM.LE.1) NSHARE(1,IJADDR) = NSHARE(1,IJADDR) + 1 NSHARE(2,IJADDR) = NSHARE(2,IJADDR) + 1 300 CONTINUE 400 CONTINUE 500 CONTINUE RETURN END