hdiff output

r32759/ido.f90 2017-06-09 15:30:14.023940838 +0100 r32758/ido.f90 2017-06-09 15:30:15.279958306 +0100
 17: !   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA 17: !   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 18: ! 18: !
 19: MODULE IDOMODULE 19: MODULE IDOMODULE
 20:      USE CONNECTDATA 20:      USE CONNECTDATA
 21:      USE CONNECTUTILS 21:      USE CONNECTUTILS
 22:      USE KEYCONNECT 22:      USE KEYCONNECT
 23:      IMPLICIT NONE 23:      IMPLICIT NONE
 24:      CONTAINS 24:      CONTAINS
 25:  25: 
 26: SUBROUTINE INITIALISE(NA,EI,Q,EF,FIN,ENDPOINTSEP) 26: SUBROUTINE INITIALISE(NA,EI,Q,EF,FIN,ENDPOINTSEP)
 27:      USE KEY, ONLY: READSP, INTIMAGE, INTLJT, INTCONSTRAINTT, FREEZENODEST, ATOMACTIVE, MLP3T, MLPVB3T,ORBITALS 27:      USE KEY, ONLY: READSP, INTIMAGE, INTLJT, INTCONSTRAINTT, FREEZENODEST, ATOMACTIVE, MLP3T, MLPVB3T
 28:      IMPLICIT NONE 28:      IMPLICIT NONE
 29:       29:      
 30:      INTEGER,INTENT(IN)           :: NA 30:      INTEGER,INTENT(IN)           :: NA
 31:      DOUBLE PRECISION,INTENT(IN)           :: ENDPOINTSEP 31:      DOUBLE PRECISION,INTENT(IN)           :: ENDPOINTSEP
 32:      DOUBLE PRECISION,POINTER              :: EI,EF 32:      DOUBLE PRECISION,POINTER              :: EI,EF
 33:      DOUBLE PRECISION,POINTER,DIMENSION(:) :: Q,FIN 33:      DOUBLE PRECISION,POINTER,DIMENSION(:) :: Q,FIN
 34:  34: 
 35:      INTEGER J2, NPLUS, NMINUS, MINPOS, NMINA, NMINB, NTSDUMP, NDUMMY, NMINDUMP, IOERROR, NMAXINT, NMININT 35:      INTEGER J2, NPLUS, NMINUS, MINPOS, NMINA, NMINB, NTSDUMP, NDUMMY, NMINDUMP, IOERROR, NMAXINT, NMININT
 36:      LOGICAL MINNEW 36:      LOGICAL MINNEW
 37:      DOUBLE PRECISION, POINTER :: ETEMP, LOCALPOINTS(:) 37:      DOUBLE PRECISION, POINTER :: ETEMP, LOCALPOINTS(:)
 39:      DOUBLE PRECISION  DJWDUMMY, CONSTRAINTE, XYZLOCAL(6*NA), GGGLOCAL(6*NA), RMSLOCAL, MINCOORDS(2,3*NA), EEELOCAL(INTIMAGE+2) 39:      DOUBLE PRECISION  DJWDUMMY, CONSTRAINTE, XYZLOCAL(6*NA), GGGLOCAL(6*NA), RMSLOCAL, MINCOORDS(2,3*NA), EEELOCAL(INTIMAGE+2)
 40:      INTEGER, ALLOCATABLE :: ACTUALMIN(:) 40:      INTEGER, ALLOCATABLE :: ACTUALMIN(:)
 41:      LOGICAL PERMUTE, IMGFREEZELOCAL(0), FREEZENODESTLOCAL 41:      LOGICAL PERMUTE, IMGFREEZELOCAL(0), FREEZENODESTLOCAL
 42: !    LOGICAL EDGEINT(INTIMAGE+1,NA,NA) 42: !    LOGICAL EDGEINT(INTIMAGE+1,NA,NA)
 43:  43: 
 44:      INTEGER INVERT,INDEX(NA),IMATCH, INTIMAGESAVE 44:      INTEGER INVERT,INDEX(NA),IMATCH, INTIMAGESAVE
 45:  45: 
 46:      NULLIFY(ETEMP, LOCALPOINTS) 46:      NULLIFY(ETEMP, LOCALPOINTS)
 47:      NATOMS=NA 47:      NATOMS=NA
 48:      NOPT=3*NATOMS 48:      NOPT=3*NATOMS
 49:      IF (MLP3T.OR.MLPVB3T.OR.ORBITALS) NOPT=NATOMS 49:      IF (MLP3T.OR.MLPVB3T) NOPT=NATOMS
 50:      ALLOCATE(G(NOPT),MI(MINRACKSIZE),TS(TSRACKSIZE)) 50:      ALLOCATE(G(NOPT),MI(MINRACKSIZE),TS(TSRACKSIZE))
 51:  51: 
 52:      ! endpoints initialisation 52:      ! endpoints initialisation
 53:      NMIN=2 53:      NMIN=2
 54:      MI(1)%DATA%E => EI 54:      MI(1)%DATA%E => EI
 55:      MI(2)%DATA%E => EF 55:      MI(2)%DATA%E => EF
 56:      MI(1)%DATA%X => Q 56:      MI(1)%DATA%X => Q
 57:      MI(2)%DATA%X => FIN 57:      MI(2)%DATA%X => FIN
 58:      ALLOCATE(MI(2)%DATA%D(1),MI(2)%DATA%NTRIES(1),MI(2)%DATA%INTNTRIES(1)) 58:      ALLOCATE(MI(2)%DATA%D(1),MI(2)%DATA%NTRIES(1),MI(2)%DATA%INTNTRIES(1))
 59:      MI(2)%DATA%D(1) = ENDPOINTSEP 59:      MI(2)%DATA%D(1) = ENDPOINTSEP


r32759/key.f90 2017-06-09 15:30:14.259944120 +0100 r32758/key.f90 2017-06-09 15:30:15.523961710 +0100
325: 325: 
326:       DOUBLE PRECISION, ALLOCATABLE ::  MLPDAT(:,:)326:       DOUBLE PRECISION, ALLOCATABLE ::  MLPDAT(:,:)
327:       INTEGER, ALLOCATABLE ::  MLPOUTCOME(:)327:       INTEGER, ALLOCATABLE ::  MLPOUTCOME(:)
328:       DOUBLE PRECISION, ALLOCATABLE ::  MLQDAT(:,:)328:       DOUBLE PRECISION, ALLOCATABLE ::  MLQDAT(:,:)
329:       INTEGER, ALLOCATABLE ::  MLQOUTCOME(:)329:       INTEGER, ALLOCATABLE ::  MLQOUTCOME(:)
330:       INTEGER, ALLOCATABLE ::  LJADDNN(:,:)330:       INTEGER, ALLOCATABLE ::  LJADDNN(:,:)
331: 331: 
332: ! OPEP variables332: ! OPEP variables
333:       LOGICAL :: OPEPT, OPEP_RNAT333:       LOGICAL :: OPEPT, OPEP_RNAT
334: 334: 
335: !Orbital variable335: 
336:       LOGICAL :: ORBITALS 
337:       INTEGER :: NROTS, NORBS, ORBVAREXPONENT 
338:       DOUBLE PRECISION, ALLOCATABLE :: R2INTS(:,:), DIPINTS(:,:,:) 
339: END MODULE KEY336: END MODULE KEY


r32759/keywords.f 2017-06-09 15:30:14.535947958 +0100 r32758/keywords.f 2017-06-09 15:30:15.787965374 +0100
 61:          USE MULTIPOT 61:          USE MULTIPOT
 62:          ! includes toggle for switching frames - required for RBAA 62:          ! includes toggle for switching frames - required for RBAA
 63:          USE MODHESS, ONLY : RBAANORMALMODET 63:          USE MODHESS, ONLY : RBAANORMALMODET
 64:          ! jdf43> for MMEINITWRAPPER 64:          ! jdf43> for MMEINITWRAPPER
 65:          USE ISO_C_BINDING, ONLY: C_NULL_CHAR 65:          USE ISO_C_BINDING, ONLY: C_NULL_CHAR
 66: ! AMBER12 66: ! AMBER12
 67:          USE AMBER12_INTERFACE_MOD, ONLY: AMBER12_ATOM, AMBER12_RESIDUE, AMBER12_ATOMS, 67:          USE AMBER12_INTERFACE_MOD, ONLY: AMBER12_ATOM, AMBER12_RESIDUE, AMBER12_ATOMS,
 68:      &                                    AMBER12_RESIDUES, AMBER12_GET_COORDS 68:      &                                    AMBER12_RESIDUES, AMBER12_GET_COORDS
 69:          USE CHIRALITY 69:          USE CHIRALITY
 70:          USE OPEP_INTERFACE_MOD, ONLY : OPEP_INIT 70:          USE OPEP_INTERFACE_MOD, ONLY : OPEP_INIT
 71:          USE ORBITALS_MOD, ONLY: ORBITALS_INIT 
 72:          use sandbox_module, only : num_atoms 71:          use sandbox_module, only : num_atoms
 73:  72: 
 74:          IMPLICIT NONE 73:          IMPLICIT NONE
 75:  74: 
 76:          DOUBLE PRECISION ::  Q(3*NATOMS) 75:          DOUBLE PRECISION ::  Q(3*NATOMS)
 77:  76: 
 78:          INTEGER NDUM, LUNIT, FUNIT, GETUNIT, MLPDATSTART, MLQDATSTART, NTYEPA 77:          INTEGER NDUM, LUNIT, FUNIT, GETUNIT, MLPDATSTART, MLQDATSTART, NTYEPA
 79:          78:         
 80:          INTEGER ITEM, NITEMS, LOC, LINE, NCR, NERROR, IR, LAST, NTYPEA, J1, J2, J3, J, I 79:          INTEGER ITEM, NITEMS, LOC, LINE, NCR, NERROR, IR, LAST, NTYPEA, J1, J2, J3, J, I
 81:          COMMON /BUFINF/ ITEM, NITEMS, LOC(132), LINE, SKIPBL, CLEAR, NCR, 80:          COMMON /BUFINF/ ITEM, NITEMS, LOC(132), LINE, SKIPBL, CLEAR, NCR,
992: 991: 
993:          MAXGAPT=.FALSE.992:          MAXGAPT=.FALSE.
994: 993: 
995:          MALONALDEHYDE=.FALSE.994:          MALONALDEHYDE=.FALSE.
996: 995: 
997:          MBPOLT= .FALSE.996:          MBPOLT= .FALSE.
998: ! OPEP stuff997: ! OPEP stuff
999:          OPEPT = .FALSE.998:          OPEPT = .FALSE.
1000:          OPEP_RNAT = .FALSE.999:          OPEP_RNAT = .FALSE.
1001: 1000: 
1002: ! cs675 ORBITALS stuff 
1003:          ORBITALS = .FALSE. 
1004:          ORBVAREXPONENT = -1 
1005:  
1006:          CLSTRINGT=.FALSE.1001:          CLSTRINGT=.FALSE.
1007:          CLSTRINGTST=.FALSE.1002:          CLSTRINGTST=.FALSE.
1008:          IF (FILTH2.EQ.0) THEN1003:          IF (FILTH2.EQ.0) THEN
1009:             OPEN (5,FILE='odata',STATUS='OLD')1004:             OPEN (5,FILE='odata',STATUS='OLD')
1010:          ELSE1005:          ELSE
1011:             WRITE(OTEMP,*) FILTH21006:             WRITE(OTEMP,*) FILTH2
1012:             WRITE(OSTRING,'(A)') 'odata.' // TRIM(ADJUSTL(OTEMP))1007:             WRITE(OSTRING,'(A)') 'odata.' // TRIM(ADJUSTL(OTEMP))
1013:             OPEN (5,FILE=OSTRING,STATUS='OLD')1008:             OPEN (5,FILE=OSTRING,STATUS='OLD')
1014:          ENDIF1009:          ENDIF
1015:          DATA_UNIT=51010:          DATA_UNIT=5
4869: ! 4864: ! 
4870:          ELSE IF (WORD.EQ.'OHCELL') THEN4865:          ELSE IF (WORD.EQ.'OHCELL') THEN
4871:             OHCELLT=.TRUE.4866:             OHCELLT=.TRUE.
4872:             BULKT=.TRUE. ! just in case we forget in odata!4867:             BULKT=.TRUE. ! just in case we forget in odata!
4873:             WRITE(*,'(A)') ' Octahedral supercell specfied'4868:             WRITE(*,'(A)') ' Octahedral supercell specfied'
4874:          ELSE IF (WORD.EQ.'ODIHE') THEN4869:          ELSE IF (WORD.EQ.'ODIHE') THEN
4875:             ODIHET=.TRUE.4870:             ODIHET=.TRUE.
4876:             WRITE(*,'(A)') 'ODIHE set: dihedral-angle order parameter will be calculated'4871:             WRITE(*,'(A)') 'ODIHE set: dihedral-angle order parameter will be calculated'
4877:             WRITE(*,'(A)') 'using the reference structure supplied in ref.crd'4872:             WRITE(*,'(A)') 'using the reference structure supplied in ref.crd'
4878: 4873: 
4879: ! Orbital localisation potential. 
4880:          ELSE IF (WORD.EQ.'ORBITALS') THEN 
4881:             ORBITALS = .TRUE. 
4882:             CALL READI(ORBVAREXPONENT) 
4883:             PRINT *, ORBVAREXPONENT 
4884:             CALL ORBITALS_INIT(NORBS, NROTS) 
4885: C4874: C
4886: C  dm368 potential4875: C  dm368 potential
4887: C4876: C
4888:       ELSE IF (WORD.EQ.'EX1D') THEN4877:       ELSE IF (WORD.EQ.'EX1D') THEN
4889:          EX1DT=.TRUE.4878:          EX1DT=.TRUE.
4890:          IF (NITEMS.GT.1) CALL READI(NATOMS)4879:          IF (NITEMS.GT.1) CALL READI(NATOMS)
4891: !         IF (MOD(NONEDAPBC,3).NE.0) THEN4880: !         IF (MOD(NONEDAPBC,3).NE.0) THEN
4892: !            WRITE(*,'(A)') 'keywords> ERROR *** lattice dimension must be a multiple of three'4881: !            WRITE(*,'(A)') 'keywords> ERROR *** lattice dimension must be a multiple of three'
4893: !            STOP4882: !            STOP
4894: !         ENDIF4883: !         ENDIF


r32759/orbitals.f90 2017-06-09 15:30:14.775951301 +0100 r32758/orbitals.f90 2017-06-09 15:30:16.015968548 +0100
  1: MODULE ORBITALS_MOD  1: svn: E195012: Unable to find repository location for 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/OPTIM/source/orbitals.f90' in revision 32758
  2:  
  3: IMPLICIT NONE 
  4:  
  5:   CONTAINS 
  6:  
  7:     SUBROUTINE ORBITALS_INIT(NORBS, NROTS) 
  8:  
  9:         ! Initialise information needed for orbital optimisation. 
 10:         ! This requires reading in the integral arrays from file. 
 11:  
 12:         USE KEY, ONLY: R2INTS, DIPINTS, ORBVAREXPONENT 
 13:  
 14:         INTEGER, INTENT(OUT) :: NORBS, NROTS 
 15:  
 16:         INTEGER :: ORB_FILE_UNIT, GETUNIT 
 17:  
 18:         CHARACTER (LEN=10) :: WORD 
 19:  
 20:         LOGICAL :: END 
 21:  
 22:         ! First check we've initialised sensibly. 
 23:         IF (ORBVAREXPONENT.LT.1) THEN 
 24:             STOP 'Need to set penalty function exponent to a positive value.' 
 25:         END IF 
 26:  
 27:         NORBS = -1 
 28:  
 29:         CALL CHECK_FILE_EXISTS('INTEGRALS_INFO', EX=.TRUE.) 
 30:  
 31:         ORB_FILE_UNIT=GETUNIT() 
 32:         OPEN(UNIT=ORB_FILE_UNIT,FILE='INTEGRALS_INFO',STATUS='OLD') 
 33:         !PRINT *,'Reading norbs...' 
 34:         READ(ORB_FILE_UNIT, *) NORBS 
 35:         NROTS = NORBS * (NORBS - 1) / 2 
 36:         !PRINT '(A,3I10)','NORBS,NROTS,ORBVAREXPONENT=',NORBS,NROTS,ORBVAREXPONENT 
 37:  
 38:         IF (NORBS.LT.1) THEN 
 39:             STOP 'Could not successfully read the system information from file.' 
 40:         ELSE IF (NROTS.LT.1) THEN 
 41:             STOP 'We appear to have no possible orbital rotations.' 
 42:         END IF 
 43:  
 44:         CLOSE(ORB_FILE_UNIT) 
 45:  
 46:         ALLOCATE(R2INTS(NORBS,NORBS)) 
 47:         ALLOCATE(DIPINTS(3,NORBS,NORBS)) 
 48:  
 49:         !PRINT '(A)','R2' 
 50:         CALL READ_INTEGRALS('INTEGRALS_R2.csv', NORBS, R2INTS(:,:)) 
 51:         !PRINT '(A)','X' 
 52:         CALL READ_INTEGRALS('INTEGRALS_X.csv', NORBS, DIPINTS(1,:,:)) 
 53:         !PRINT '(A)','Y' 
 54:         CALL READ_INTEGRALS('INTEGRALS_Y.csv', NORBS, DIPINTS(2,:,:)) 
 55:         !PRINT '(A)','Z' 
 56:         CALL READ_INTEGRALS('INTEGRALS_Z.csv', NORBS, DIPINTS(3,:,:)) 
 57:  
 58:     END SUBROUTINE ORBITALS_INIT 
 59:  
 60:     SUBROUTINE READ_INTEGRALS(INT_FILENAME, NORBS, INTS) 
 61:  
 62:         ! Given filename of intergral files reads the integrals within 
 63:         ! it into the array INTS. Integrals are assumed to always have 
 64:         ! dimension NORBSxNORBS. 
 65:  
 66:         CHARACTER (LEN=*), INTENT(IN) :: INT_FILENAME 
 67:         INTEGER, INTENT(IN) :: NORBS 
 68:         DOUBLE PRECISION, INTENT(OUT) :: INTS(NORBS, NORBS) 
 69:         CHARACTER (LEN=32) :: CALCID 
 70:         INTEGER :: INT_FILE_UNIT, I, GETUNIT 
 71:         LOGICAL :: END 
 72:  
 73:         CALL CHECK_FILE_EXISTS(INT_FILENAME, EX=.TRUE.) 
 74:         INT_FILE_UNIT=GETUNIT() 
 75:         OPEN(UNIT=INT_FILE_UNIT,FILE=INT_FILENAME,STATUS='OLD') 
 76: !        CALL FILE_OPEN(INT_FILENAME,INT_FILE_UNIT,.FALSE.) 
 77: !        CALL INPUT(END, INT_FILE_UNIT) 
 78:         READ(INT_FILE_UNIT, '(A)') CALCID 
 79:         !PRINT '(A,I10,X,A)','INT_FILE_UNIT,CALCID=',INT_FILE_UNIT,CALCID 
 80:         !PRINT '(A,I10,A)','NORBS=',NORBS 
 81:  
 82:         DO I = 1, NORBS 
 83:            READ(INT_FILE_UNIT,*) INTS(I,1:NORBS)  
 84:            !PRINT *,INTS(I,1:NORBS) 
 85:         END DO 
 86:  
 87:         CLOSE(INT_FILE_UNIT) 
 88:  
 89:     END SUBROUTINE READ_INTEGRALS 
 90:  
 91:     SUBROUTINE GET_ORBITAL_LOCALITY(X, GRAD, LOCALITY, GTEST, SECT) 
 92:  
 93:         ! Obtains the Power of the Orbital Variance locality measure for 
 94:         ! a given position. If GTEST is set also returns the gradient. 
 95:  
 96:         ! Note that for intermediates dependent upon coefficient of current 
 97:         ! (rotated) orbital set in the original MO basis we use indexing 
 98:         ! of (ORIG_ORB,NEW_ORB) in all cases. 
 99:         ! For 3-index quantities we use (ORIG_ORB1,ORIG_ORB2,NEW_ORB) 
100:  
101:         USE KEY, ONLY: NROTS, NORBS, R2INTS, DIPINTS, ORBVAREXPONENT 
102:         !USE MODHESS 
103:  
104:         DOUBLE PRECISION, INTENT(IN) :: X(NROTS) 
105:         DOUBLE PRECISION, INTENT(OUT) :: GRAD(NROTS), LOCALITY 
106:         LOGICAL, INTENT(IN) :: GTEST, SECT 
107:  
108:         DOUBLE PRECISION :: ROTATION(NORBS,NORBS), ORBVAR(NORBS) 
109:         DOUBLE PRECISION :: ORBDIPOLES(3,NORBS), IND_ORBVAR_DERIV(NORBS,NORBS) 
110:         DOUBLE PRECISION :: PENALTY_DERIV(NORBS,NORBS), ROT_DERIV(NORBS,NORBS) 
111:         DOUBLE PRECISION :: PROD_DERIVS(NORBS,NORBS) 
112:  
113:         INTEGER :: I,J,K 
114:  
115:         CALL GET_ROTATION(X, ROTATION) 
116:  
117:         CALL GET_CURRENT_ORB_DIPOLES(ROTATION, ORBDIPOLES) 
118:         DO I = 1, NORBS 
119:             ORBVAR(I) = -SUM(ORBDIPOLES(1:3,I)**2) 
120:             DO J = 1, NORBS 
121:                 DO K = 1, NORBS 
122:                     ORBVAR(I)=ORBVAR(I)+R2INTS(J,K)*ROTATION(J,I)*ROTATION(K,I) 
123:                 END DO 
124:             END DO 
125:             !PRINT *,'ORBITAL',I,',ORBVAR=',ORBVAR(I),',DIPOLES=',ORBDIPOLES(1:3,I) 
126:         END DO 
127:  
128:         !PRINT *,'ORBVARS=',ORBVAR 
129:         !PRINT *,'ORBVARSEXP=',ORBVAR**ORBVAREXPONENT 
130:  
131:         LOCALITY = SUM(ORBVAR ** ORBVAREXPONENT) 
132:  
133:         !PRINT *,'LOCALITY=',LOCALITY 
134:  
135:         IF (GTEST) THEN 
136:             !PRINT *,'ORBDIPOLES=' 
137:             !DO J = 1,3 
138:             !    PRINT *,ORBDIPOLES(J,1:NORBS) 
139:             !END DO 
140:             CALL GET_INDIVIDUAL_DERIVATIVES(ROTATION, ORBDIPOLES, IND_ORBVAR_DERIV) 
141:             !PRINT *,'IND_ORBVAR_DERIV=' 
142:             !DO J = 1,NORBS 
143:             !    PRINT *, IND_ORBVAR_DERIV(J,1:NORBS) 
144:             !END DO 
145:  
146:             CALL GET_PENALTY_DERIVATIVES(IND_ORBVAR_DERIV, ORBVAR, PENALTY_DERIV) 
147:             !PRINT *,'PEN_DERIV=' 
148:             !DO J = 1,NORBS 
149:             !    PRINT *, PENALTY_DERIV(J,1:NORBS) 
150:             !END DO 
151:  
152:             DO I = 1, NROTS 
153:                 CALL GET_ROTATION_DERIVATIVE(X, I, ROT_DERIV) 
154:                 !PRINT *,'ROT_DERIV',I,'=' 
155:                 !DO J = 1,NORBS 
156:                 !    PRINT *, ROT_DERIV(J,1:NORBS) 
157:                 !END DO 
158:  
159:                 CALL ELEMENTWISE_MULTIPLICATION(ROT_DERIV,PENALTY_DERIV,PROD_DERIVS) 
160:                 !PRINT *,'OVERALL_PROD=' 
161:                 !DO J = 1,NORBS 
162:                 !    PRINT *, PROD_DERIVS(J,:) 
163:                 !END DO 
164:                 GRAD(I) = SUM(PROD_DERIVS) 
165:                 !PRINT *,'GRAD',I,'=',GRAD(I) 
166:             END DO 
167:  
168:         END IF 
169:  
170:         IF (SECT) THEN 
171:             ! If haven't calculated gradient as well need to calculate additional intermediates. 
172:             IF (.NOT.GTEST) THEN 
173:                 CALL GET_INDIVIDUAL_DERIVATIVES(ROTATION, ORBDIPOLES, IND_ORBVAR_DERIV) 
174:                 CALL GET_PENALTY_DERIVATIVES(IND_ORBVAR_DERIV, ORBVAR, PENALTY_DERIV) 
175:             END IF 
176:  
177:             !IF (.NOT.ALLOCATED(HESS)) THEN 
178:             !    PRINT *, 'Allocating Hessian with size ',NROTS,'x',NROTS 
179:             !    ALLOCATE(HESS(NROTS,NROTS)) 
180:             !END IF 
181:             !PRINT*,'hessian bounds:',LBOUND(HESS,DIM=1),UBOUND(HESS,DIM=1),LBOUND(HESS,DIM=2),UBOUND(HESS,DIM=2) 
182:             !PRINT*,'hessian size:',SIZE(HESS),SHAPE(HESS) 
183:             CALL CALC_HESSIAN(X,ROTATION,ORBVAR,ORBDIPOLES,IND_ORBVAR_DERIV,PENALTY_DERIV) 
184:  
185:         END IF 
186:  
187:     END SUBROUTINE GET_ORBITAL_LOCALITY 
188:  
189:     SUBROUTINE ELEMENTWISE_MULTIPLICATION(MATA,MATB,MATC) 
190:  
191:         USE KEY, ONLY: NORBS 
192:  
193:         DOUBLE PRECISION, INTENT(IN) :: MATA(NORBS,NORBS), MATB(NORBS,NORBS) 
194:  
195:         DOUBLE PRECISION, INTENT(OUT) :: MATC(NORBS, NORBS) 
196:  
197:         INTEGER :: I, J 
198:  
199:         DO I = 1, NORBS 
200:             DO J = 1, NORBS 
201:                 MATC(I,J) = MATA(I,J) * MATB(I,J) 
202:             END DO 
203:         END DO 
204:  
205:     END SUBROUTINE ELEMENTWISE_MULTIPLICATION 
206:  
207:     SUBROUTINE GET_ROTATION(COORDS, ROTATION) 
208:  
209:         ! Obtains the overall rotation resulting from the direct product 
210:         ! of the Givens rotations corresponding to the specified coordinates. 
211:  
212:         USE KEY, ONLY: NROTS, NORBS 
213:  
214:         DOUBLE PRECISION, INTENT(IN) :: COORDS(NROTS) 
215:  
216:         DOUBLE PRECISION, INTENT(OUT) :: ROTATION(NORBS,NORBS) 
217:         DOUBLE PRECISION :: NEXT_ROT(NORBS,NORBS) 
218:  
219:         INTEGER :: I, J, TRIIND 
220:  
221:         ROTATION(1:NORBS,1:NORBS) = 0.0 
222:  
223:         DO TRIIND = 1, NORBS 
224:             ROTATION(TRIIND, TRIIND) = 1.0 
225:         END DO 
226:  
227:         DO TRIIND = 1, NROTS 
228:             CALL DECODE_TRIIND(TRIIND, I, J) 
229:             CALL GET_GIVENS_ROTATION(I,J,COORDS(TRIIND),NEXT_ROT) 
230:             ROTATION = MATMUL(ROTATION,NEXT_ROT) 
231:         END DO 
232:  
233:     END SUBROUTINE GET_ROTATION 
234:  
235:     SUBROUTINE GET_ROTATION_DERIVATIVE(COORDS, DERIVTRIIND, DERIV_ROTATION) 
236:  
237:         ! Obtain the derivative of the rotation matrix with respect to a single 
238:         ! rotation angle theta at the current coordinates. 
239:  
240:         USE KEY, ONLY: NROTS, NORBS 
241:  
242:         DOUBLE PRECISION, INTENT(IN) :: COORDS(NROTS) 
243:         INTEGER, INTENT(IN) :: DERIVTRIIND 
244:  
245:         DOUBLE PRECISION, INTENT(OUT) :: DERIV_ROTATION(NORBS,NORBS) 
246:  
247:         DOUBLE PRECISION :: NEXT_ROT(NORBS,NORBS) 
248:         INTEGER :: TRIIND, I, J 
249:  
250:         DERIV_ROTATION(1:NORBS,1:NORBS) = 0.0 
251:  
252:         DO TRIIND = 1, NORBS 
253:             DERIV_ROTATION(TRIIND, TRIIND) = 1.0 
254:         END DO 
255:  
256:         DO TRIIND = 1, NROTS 
257:             CALL DECODE_TRIIND(TRIIND, I, J) 
258:             IF (TRIIND.EQ.DERIVTRIIND) THEN 
259:                 CALL GET_DERIVATIVE_GIVENS_ROTATION(I,J,COORDS(TRIIND),NEXT_ROT) 
260:             ELSE 
261:                 CALL GET_GIVENS_ROTATION(I,J,COORDS(TRIIND),NEXT_ROT) 
262:             END IF 
263:             DERIV_ROTATION = MATMUL(DERIV_ROTATION,NEXT_ROT) 
264:         END DO 
265:  
266:         !DO TRIIND = 1, NORBS 
267:         !    PRINT *,TRIIND,'DERIV_ROT=',DERIV_ROTATION(TRIIND,1:NORBS) 
268:         !END DO 
269:  
270:     END SUBROUTINE GET_ROTATION_DERIVATIVE 
271:  
272:     SUBROUTINE GET_ROTATION_SECOND_DERIVATIVE(COORDS, DERIVTRIINDA, DERIVTRIINDB, SECDERIV_ROTATION) 
273:  
274:         ! Obtain the derivative of the rotation matrix with respect to a single 
275:         ! rotation angle theta at the current coordinates. 
276:  
277:         USE KEY, ONLY: NROTS, NORBS 
278:  
279:         DOUBLE PRECISION, INTENT(IN) :: COORDS(NROTS) 
280:         INTEGER, INTENT(IN) :: DERIVTRIINDA, DERIVTRIINDB 
281:  
282:         DOUBLE PRECISION, INTENT(OUT) :: SECDERIV_ROTATION(NORBS,NORBS) 
283:  
284:         DOUBLE PRECISION :: NEXT_ROT(NORBS,NORBS) 
285:         INTEGER :: TRIIND, I, J 
286:  
287:         SECDERIV_ROTATION(1:NORBS,1:NORBS) = 0.0 
288:  
289:         DO TRIIND = 1, NORBS 
290:             SECDERIV_ROTATION(TRIIND, TRIIND) = 1.0 
291:         END DO 
292:  
293:         DO TRIIND = 1, NROTS 
294:             CALL DECODE_TRIIND(TRIIND, I, J) 
295:             IF (TRIIND.EQ.DERIVTRIINDA) THEN 
296:                 IF (DERIVTRIINDA.EQ.DERIVTRIINDB) THEN 
297:                     CALL GET_SECOND_DERIVATIVE_GIVENS_ROTATION(I,J,COORDS(TRIIND),NEXT_ROT) 
298:                 ELSE 
299:                     CALL GET_DERIVATIVE_GIVENS_ROTATION(I,J,COORDS(TRIIND),NEXT_ROT) 
300:                 END IF 
301:             ELSE IF (TRIIND.EQ.DERIVTRIINDB) THEN 
302:                 CALL GET_DERIVATIVE_GIVENS_ROTATION(I,J,COORDS(TRIIND),NEXT_ROT) 
303:             ELSE 
304:                 CALL GET_GIVENS_ROTATION(I,J,COORDS(TRIIND),NEXT_ROT) 
305:             END IF 
306:             SECDERIV_ROTATION = MATMUL(SECDERIV_ROTATION,NEXT_ROT) 
307:         END DO 
308:  
309:         !DO TRIIND = 1, NORBS 
310:         !    PRINT *,TRIIND,'DERIV_ROT=',SECDERIV_ROTATION(TRIIND,1:NORBS) 
311:         !END DO 
312:  
313:     END SUBROUTINE GET_ROTATION_SECOND_DERIVATIVE 
314:  
315:     SUBROUTINE GET_GIVENS_ROTATION(I, J, THETA, MAT) 
316:  
317:         ! Obtains Givens rotation between orbitals I and J with rotation angle theta. 
318:  
319:         USE KEY, ONLY: NORBS 
320:  
321:         INTEGER, INTENT(IN) :: I, J 
322:         DOUBLE PRECISION, INTENT(IN) :: THETA 
323:         DOUBLE PRECISION, INTENT(OUT) :: MAT(NORBS,NORBS) 
324:  
325:         INTEGER :: VAL 
326:  
327:         MAT(1:NORBS,1:NORBS) = 0.0 
328:  
329:         DO VAL = 1, NORBS 
330:             MAT(VAL, VAL) = 1.0 
331:         END DO 
332:  
333:         MAT(I,I) = COS(THETA) 
334:         MAT(J,J) = COS(THETA) 
335:         MAT(I,J) = -SIN(THETA) 
336:         MAT(J,I) = SIN(THETA) 
337:  
338:         !PRINT *, 'GIVENSROT',I,'-',J,'=' 
339:         !DO VAL=1,NORBS 
340:         !    PRINT *,MAT(VAL,1:NORBS) 
341:         !END DO 
342:  
343:     END SUBROUTINE GET_GIVENS_ROTATION 
344:  
345:     SUBROUTINE GET_DERIVATIVE_GIVENS_ROTATION(I, J, THETA, MAT) 
346:  
347:         ! Obtains the first derivative of the Givens rotation between orbitals I and J 
348:         !  wrt the rotation angle at angle theta as a matrix. 
349:  
350:         USE KEY, ONLY: NORBS 
351:  
352:         INTEGER, INTENT(IN) :: I, J 
353:         DOUBLE PRECISION, INTENT(IN) :: THETA 
354:         DOUBLE PRECISION, INTENT(OUT) :: MAT(NORBS,NORBS) 
355:  
356:         INTEGER :: VAL 
357:  
358:         MAT(1:NORBS,1:NORBS) = 0.0 
359:  
360:         MAT(I,I) = -SIN(THETA) 
361:         MAT(J,J) = -SIN(THETA) 
362:         MAT(I,J) = -COS(THETA) 
363:         MAT(J,I) = COS(THETA) 
364:         !PRINT *, 'DERIVGIVENSROT',I,'-',J,'=' 
365:         !DO VAL=1,NORBS 
366:         !    PRINT *,MAT(VAL,1:NORBS) 
367:         !END DO 
368:  
369:     END SUBROUTINE GET_DERIVATIVE_GIVENS_ROTATION 
370:  
371:     SUBROUTINE GET_SECOND_DERIVATIVE_GIVENS_ROTATION(I, J, THETA, MAT) 
372:  
373:         ! Obtains the second derivative of the Givens rotation between orbitals I and J 
374:         !  wrt the rotation angle at angle theta as a matrix. 
375:  
376:         USE KEY, ONLY: NORBS 
377:  
378:         INTEGER, INTENT(IN) :: I, J 
379:         DOUBLE PRECISION, INTENT(IN) :: THETA 
380:         DOUBLE PRECISION, INTENT(OUT) :: MAT(NORBS,NORBS) 
381:  
382:         INTEGER :: VAL 
383:  
384:         MAT(1:NORBS,1:NORBS) = 0.0 
385:  
386:         MAT(I,I) = -COS(THETA) 
387:         MAT(J,J) = -COS(THETA) 
388:         MAT(I,J) = SIN(THETA) 
389:         MAT(J,I) = -SIN(THETA) 
390:         !PRINT *, 'SECONDDERIVGIVENSROT',I,'-',J,'=' 
391:         !DO VAL=1,NORBS 
392:         !    PRINT *,MAT(VAL,1:NORBS) 
393:         !END DO 
394:  
395:     END SUBROUTINE GET_SECOND_DERIVATIVE_GIVENS_ROTATION 
396:  
397:     SUBROUTINE GET_CURRENT_ORB_DIPOLES(ROTATION, ORBDIPOLES) 
398:  
399:         ! Obtain the dipole moments of the orbital set created by the current rotation of the starting MOs. 
400:  
401:         USE KEY, ONLY: NORBS, DIPINTS 
402:  
403:         DOUBLE PRECISION, INTENT(IN) :: ROTATION(NORBS, NORBS) 
404:         DOUBLE PRECISION :: ORBDIPOLES(3,NORBS) 
405:  
406:         INTEGER :: CURRENT_ORB, ORIG_ORBA, ORIG_ORBB, DIRECTION 
407:  
408:         ORBDIPOLES(1:3,1:NORBS) = 0.0 
409:  
410:         DO DIRECTION = 1,3 
411:             DO CURRENT_ORB = 1, NORBS 
412:                 DO ORIG_ORBA = 1, NORBS 
413:                     DO ORIG_ORBB = 1, NORBS 
414:                         ORBDIPOLES(DIRECTION,CURRENT_ORB) = ORBDIPOLES(DIRECTION,CURRENT_ORB) + & 
415:                             ROTATION(ORIG_ORBA,CURRENT_ORB) * ROTATION(ORIG_ORBB,CURRENT_ORB) * & 
416:                             DIPINTS(DIRECTION,ORIG_ORBA,ORIG_ORBB) 
417:                     END DO 
418:                 END DO 
419:             END DO 
420:         END DO 
421:  
422:         ! Comparison using matrix multiplication. 
423:         !DO DIRECTION = 1,3 
424:         !    PRINT *, 'Prev and matmul dipoles, direction', DIRECTION, '.' 
425:         !    TEMPCURRENTDIPOLES = MATMUL(MATMUL(TRANSPOSE(ROTATION),DIPINTS(DIRECTION,1:NORBS,1:NORBS)),ROTATION) 
426:         !    DO CURRENT_ORB = 1, NORBS 
427:         !        PRINT *, ORBDIPOLES(DIRECTION, CURRENT_ORB), TEMPCURRENTDIPOLES(CURRENT_ORB, CURRENT_ORB) 
428:         !    END DO 
429:         !END DO 
430:  
431:     END SUBROUTINE GET_CURRENT_ORB_DIPOLES 
432:  
433:     SUBROUTINE DECODE_TRIIND(TRIIND, I, J) 
434:  
435:         ! Given a triangular index obtains the 2D indices corresponding to it. 
436:  
437:         USE KEY, ONLY: NORBS 
438:  
439:         INTEGER, INTENT(IN) :: TRIIND 
440:         INTEGER, INTENT(OUT) :: I, J 
441:  
442:         INTEGER :: RUNNING_SUM, VAL, I_, J_ 
443:  
444:         RUNNING_SUM = 0 
445:  
446:         DO VAL = 1, NORBS 
447:             IF (RUNNING_SUM + VAL .GE. TRIIND) THEN 
448:                 I_ = VAL + 1 
449:                 J_ = TRIIND - RUNNING_SUM 
450:                 EXIT 
451:             ELSE 
452:                 RUNNING_SUM = RUNNING_SUM + VAL 
453:             END IF 
454:  
455:         END DO 
456:  
457:         IF (I_ < J_) THEN 
458:             I = J_ 
459:             J = I_ 
460:         ELSE 
461:             I = I_ 
462:             J = J_ 
463:         END IF 
464:  
465:         !PRINT '(A,3I10)','TRIIND,I,J=',TRIIND,I,J 
466:  
467:     END SUBROUTINE DECODE_TRIIND 
468:  
469:     SUBROUTINE GET_INDIVIDUAL_DERIVATIVES(CURRENT_ROTATION, CURRENT_DIP, INDIVIDUAL_DERIVATIVES) 
470:  
471:         ! Generates the matrix of all derivatives of different orbital variances w.r.t. 
472:         ! the coefficients of their expansion in the original MO basis. 
473:  
474:         USE KEY, ONLY: NORBS 
475:  
476:         DOUBLE PRECISION, INTENT(IN) :: CURRENT_ROTATION(NORBS, NORBS), CURRENT_DIP(3,NORBS) 
477:  
478:         DOUBLE PRECISION, INTENT(OUT) :: INDIVIDUAL_DERIVATIVES(NORBS,NORBS) 
479:  
480:         INTEGER :: ORIGORB, NEWORB 
481:  
482:         !PRINT *,'CURRENTROTATION=' 
483:         !DO NEWORB=1,NORBS 
484:         !    PRINT*,CURRENT_ROTATION(NEWORB,1:NORBS) 
485:         !END DO 
486:  
487:         DO NEWORB = 1, NORBS 
488:             DO ORIGORB = 1, NORBS 
489:                 CALL GET_SINGLE_ORBVAR_DERIVATIVE(CURRENT_ROTATION, CURRENT_DIP, NEWORB, ORIGORB, INDIVIDUAL_DERIVATIVES(ORIGORB,NEWORB)) 
490:             END DO 
491:         END DO 
492:  
493:     END SUBROUTINE GET_INDIVIDUAL_DERIVATIVES 
494:  
495:     SUBROUTINE GET_SINGLE_ORBVAR_DERIVATIVE(CURRENT_ROTATION, CURRENT_DIP, NEW_ORB, ORIG_ORB, DERIVATIVE) 
496:  
497:         ! Obtains the derivative of a single orbital at the current rotation's orbital 
498:         ! variance w.r.t. to the coefficient an original (unrotated) orbital. 
499:  
500:         USE KEY, ONLY: NORBS, R2INTS, DIPINTS 
501:  
502:         DOUBLE PRECISION, INTENT(IN) :: CURRENT_ROTATION(NORBS, NORBS), CURRENT_DIP(3,NORBS) 
503:  
504:         DOUBLE PRECISION, INTENT(OUT) :: DERIVATIVE 
505:  
506:         INTEGER,INTENT(IN) :: NEW_ORB, ORIG_ORB 
507:  
508:         INTEGER :: OTHER_ORIG_ORB 
509:  
510:         DOUBLE PRECISION :: CONTRIB 
511:  
512:         DERIVATIVE = 0.0 
513:  
514:         DO OTHER_ORIG_ORB = 1, NORBS 
515:             CONTRIB = 2 * CURRENT_ROTATION(OTHER_ORIG_ORB,NEW_ORB)*R2INTS(ORIG_ORB, OTHER_ORIG_ORB) 
516:  
517:             DERIVATIVE = DERIVATIVE + CONTRIB 
518:  
519:  
520:             CONTRIB = 4 * CURRENT_ROTATION(OTHER_ORIG_ORB,NEW_ORB)*& 
521:                         (CURRENT_DIP(1,NEW_ORB)*DIPINTS(1,ORIG_ORB,OTHER_ORIG_ORB)& 
522:                         +CURRENT_DIP(2,NEW_ORB)*DIPINTS(2,ORIG_ORB,OTHER_ORIG_ORB)& 
523:                         +CURRENT_DIP(3,NEW_ORB)*DIPINTS(3,ORIG_ORB,OTHER_ORIG_ORB)) 
524:  
525:             DERIVATIVE = DERIVATIVE - CONTRIB 
526:  
527:         END DO 
528:  
529:     END SUBROUTINE GET_SINGLE_ORBVAR_DERIVATIVE 
530:  
531:     SUBROUTINE GET_PENALTY_DERIVATIVES(IND_ORBVAR_DERIV, ORB_VAR, PENALTY_DERIVATIVES) 
532:  
533:         ! Get the derivative of the overall penalty function with respect to each element 
534:         ! of the rotation matrix. 
535:  
536:         USE KEY, ONLY: NORBS, ORBVAREXPONENT 
537:  
538:         DOUBLE PRECISION, INTENT(IN) :: IND_ORBVAR_DERIV(NORBS,NORBS), ORB_VAR(NORBS) 
539:         DOUBLE PRECISION, INTENT(OUT) :: PENALTY_DERIVATIVES(NORBS, NORBS) 
540:  
541:         INTEGER :: ORIG_ORB, NEW_ORB 
542:  
543:         PENALTY_DERIVATIVES = 0.0 
544:  
545:         DO ORIG_ORB = 1, NORBS 
546:             DO NEW_ORB = 1, NORBS 
547:                 PENALTY_DERIVATIVES(ORIG_ORB,NEW_ORB) = & 
548:                     ORBVAREXPONENT * IND_ORBVAR_DERIV(ORIG_ORB,NEW_ORB) * ORB_VAR(NEW_ORB) ** (ORBVAREXPONENT-1) 
549:             END DO 
550:         END DO 
551:  
552:     END SUBROUTINE GET_PENALTY_DERIVATIVES 
553:  
554:     SUBROUTINE CALC_HESSIAN(COORDS,ROTATION,ORBVAR,ORBDIPOLES,IND_ORBVAR_DERIV,PENALTY_DERIV) 
555:  
556:         USE KEY, ONLY: NROTS, NORBS 
557:         USE MODHESS 
558:  
559:         DOUBLE PRECISION, INTENT(IN) :: COORDS(NROTS) 
560:         DOUBLE PRECISION, INTENT(IN) :: ROTATION(NORBS,NORBS), ORBDIPOLES(3,NORBS), IND_ORBVAR_DERIV(NORBS,NORBS) 
561:         DOUBLE PRECISION, INTENT(IN) :: ORBVAR(NORBS), PENALTY_DERIV(NORBS,NORBS) 
562:  
563:         DOUBLE PRECISION :: IND_ORBVAR_SECOND_DERIV(NORBS,NORBS,NORBS) 
564:         DOUBLE PRECISION :: PENALTY_SECOND_DERIV(NORBS,NORBS,NORBS), TEMP 
565:         DOUBLE PRECISION :: ROTMATA(NORBS,NORBS), ROTMATB(NORBS,NORBS) 
566:  
567:         INTEGER :: TRIINDA, TRIINDB, NEW_ORB, ORIG_ORBA, ORIG_ORBB 
568:  
569:         ! Gradually work our way up the chain rule... 
570:         CALL CALC_ORBVAR_SECOND_DERIV(ROTATION,ORBDIPOLES,IND_ORBVAR_SECOND_DERIV) 
571:  
572:         CALL CALC_PENALTY_SECOND_DERIV(ORBVAR,IND_ORBVAR_DERIV,IND_ORBVAR_SECOND_DERIV,PENALTY_SECOND_DERIV) 
573:  
574:         !PRINT*,'hessian bounds:',LBOUND(HESS,DIM=1),UBOUND(HESS,DIM=1),LBOUND(HESS,DIM=2),UBOUND(HESS,DIM=2) 
575:         !PRINT*,'hessian size:',SIZE(HESS),SHAPE(HESS) 
576:  
577:         !PRINT*,NROTS,NORBS 
578:  
579:         !PRINT*,'IND_ORBVAR_SECOND_DERIV=' 
580:         !DO TRIINDA=1,NORBS 
581:         !    DO TRIINDB=1,NORBS 
582:         !        PRINT*,TRIINDA,TRIINDB,IND_ORBVAR_SECOND_DERIV(TRIINDA,TRIINDB,1:NORBS) 
583:         !    END DO 
584:         !END DO 
585:         !PRINT*,'PENALTY_SECOND_DERIV=' 
586:         !DO TRIINDA=1,NORBS 
587:         !    DO TRIINDB=1,NORBS 
588:         !        PRINT*,TRIINDA,TRIINDB,PENALTY_SECOND_DERIV(TRIINDA,TRIINDB,1:NORBS) 
589:         !    END DO 
590:         !END DO 
591:  
592:         HESS(1:NROTS,1:NROTS) = 0.0 
593:  
594:         ! Now calculate actual hessian! 
595:         DO TRIINDA=1,NROTS 
596:             DO TRIINDB=1,NROTS 
597:                 ! First deal with term proportional to penalty first deriv wrt coeffs 
598:                 CALL GET_ROTATION_SECOND_DERIVATIVE(COORDS,TRIINDA,TRIINDB,ROTMATA) 
599:  
600:                 ! Use ROTMATB as scratch space. 
601:                 CALL ELEMENTWISE_MULTIPLICATION(ROTMATA(1:NORBS,1:NORBS),PENALTY_DERIV(1:NORBS,1:NORBS),ROTMATB) 
602:                 !PRINT*,'PRODMAT SECDERIV=' 
603:                 !DO NEW_ORB = 1,NORBS 
604:                 !    PRINT*,ROTMATB(NEW_ORB,1:NORBS) 
605:                 !END DO 
606:  
607:                 HESS(TRIINDA,TRIINDB) = SUM(ROTMATB) 
608:                 ! Now deal with term proportional to penalty second deriv wrt coeffs. 
609:                 CALL GET_ROTATION_DERIVATIVE(COORDS,TRIINDA,ROTMATA) 
610:                 CALL GET_ROTATION_DERIVATIVE(COORDS,TRIINDB,ROTMATB) 
611:                 DO NEW_ORB = 1,NORBS 
612:                     DO ORIG_ORBA = 1,NORBS 
613:                         DO ORIG_ORBB = 1,NORBS 
614:                             HESS(TRIINDA,TRIINDB) = HESS(TRIINDA,TRIINDB) + PENALTY_SECOND_DERIV(ORIG_ORBA,ORIG_ORBB,NEW_ORB)*ROTMATA(ORIG_ORBA,NEW_ORB)*ROTMATB(ORIG_ORBB,NEW_ORB) 
615:                         END DO 
616:                     END DO 
617:                 END DO 
618:             END DO 
619:         END DO 
620:          
621:  
622:     END SUBROUTINE CALC_HESSIAN 
623:  
624:     SUBROUTINE CALC_ORBVAR_SECOND_DERIV(ROTATION,ORBDIPOLES,IND_ORBVAR_SECOND_DERIV) 
625:  
626:         USE KEY, ONLY: NORBS, R2INTS, DIPINTS 
627:  
628:         DOUBLE PRECISION, INTENT(IN) :: ROTATION(NORBS,NORBS), ORBDIPOLES(3,NORBS) 
629:         DOUBLE PRECISION, INTENT(OUT) :: IND_ORBVAR_SECOND_DERIV(NORBS,NORBS,NORBS) 
630:  
631:         DOUBLE PRECISION :: INTERMEDIATE_SUM(NORBS) 
632:  
633:         INTEGER :: ORIG_ORBA, ORIG_ORBB, NEW_ORB, DIRECTION 
634:  
635:         ! First set values according to [x2] term. 
636:         DO ORIG_ORBA = 1,NORBS 
637:             DO ORIG_ORBB = 1,NORBS 
638:                 IND_ORBVAR_SECOND_DERIV(ORIG_ORBA,ORIG_ORBB,1:NORBS) = 2 * R2INTS(ORIG_ORBA,ORIG_ORBB) 
639:             END DO 
640:         END DO 
641:  
642:         DO NEW_ORB = 1, NORBS 
643:             DO DIRECTION = 1,3 
644:                 ! Calc \sum_i C_ip [x]_ij as for each j value (and given p). 
645:                 INTERMEDIATE_SUM(1:NORBS) = 0.0 
646:                 DO ORIG_ORBA = 1,NORBS 
647:                     DO ORIG_ORBB = 1,NORBS 
648:                         INTERMEDIATE_SUM(ORIG_ORBA) = INTERMEDIATE_SUM(ORIG_ORBA) + ROTATION(ORIG_ORBB,NEW_ORB) * DIPINTS(DIRECTION,ORIG_ORBA,ORIG_ORBB) 
649:                     END DO 
650:                 END DO 
651:                 !PRINT*,'DIRECTION=',DIRECTION 
652:                 !PRINT*,'DIPINTS=',DIPINTS(DIRECTION,:,:) 
653:                 !PRINT*,'INTERMEDIATE_SUM=',INTERMEDIATE_SUM 
654:                 DO ORIG_ORBA = 1,NORBS 
655:                     DO ORIG_ORBB = 1,NORBS 
656:                         IND_ORBVAR_SECOND_DERIV(ORIG_ORBA,ORIG_ORBB,NEW_ORB) = IND_ORBVAR_SECOND_DERIV(ORIG_ORBA,ORIG_ORBB,NEW_ORB) & 
657:                                             - 4*INTERMEDIATE_SUM(ORIG_ORBA)*INTERMEDIATE_SUM(ORIG_ORBB) & 
658:                                             - 2*ORBDIPOLES(DIRECTION,NEW_ORB)*DIPINTS(DIRECTION,ORIG_ORBA,ORIG_ORBB) 
659:                     END DO 
660:                 END DO 
661:  
662:             END DO 
663:         END DO 
664:  
665:     END SUBROUTINE CALC_ORBVAR_SECOND_DERIV 
666:  
667:     SUBROUTINE CALC_PENALTY_SECOND_DERIV(ORBVAR,IND_ORBVAR_DERIV,IND_ORBVAR_SECOND_DERIV,PENALTY_SECOND_DERIV) 
668:  
669:         USE KEY, ONLY: ORBVAREXPONENT, NORBS 
670:  
671:         DOUBLE PRECISION, INTENT(IN) :: ORBVAR(NORBS), IND_ORBVAR_DERIV(NORBS,NORBS) 
672:  
673:         DOUBLE PRECISION, INTENT(IN) :: IND_ORBVAR_SECOND_DERIV(NORBS,NORBS,NORBS) 
674:         DOUBLE PRECISION, INTENT(OUT) :: PENALTY_SECOND_DERIV(NORBS,NORBS,NORBS) 
675:  
676:         INTEGER :: NEW_ORB, ORIG_ORBA, ORIG_ORBB 
677:  
678:         DO NEW_ORB = 1,NORBS 
679:             DO ORIG_ORBA = 1,NORBS 
680:                 DO ORIG_ORBB = 1,NORBS 
681:                     PENALTY_SECOND_DERIV(ORIG_ORBA,ORIG_ORBB,NEW_ORB) = & 
682:                                 ORBVAREXPONENT * IND_ORBVAR_SECOND_DERIV(ORIG_ORBA,ORIG_ORBB,NEW_ORB)*(ORBVAR(NEW_ORB))**(ORBVAREXPONENT-1) 
683:                     IF (ORBVAREXPONENT.GT.1) THEN 
684:                         PENALTY_SECOND_DERIV(ORIG_ORBA,ORIG_ORBB,NEW_ORB) = PENALTY_SECOND_DERIV(ORIG_ORBA,ORIG_ORBB,NEW_ORB)& 
685:                                 +ORBVAREXPONENT*(ORBVAREXPONENT-1)*IND_ORBVAR_DERIV(ORIG_ORBA,NEW_ORB)*IND_ORBVAR_DERIV(ORIG_ORBB,NEW_ORB)*(ORBVAR(NEW_ORB))**(ORBVAREXPONENT-2) 
686:                     END IF 
687:                 END DO 
688:             END DO 
689:         END DO 
690:  
691:     END SUBROUTINE CALC_PENALTY_SECOND_DERIV 
692:  
693:     SUBROUTINE ORBITALS_FINISH() 
694:  
695:         ! Deallocate arrays in commons. 
696:  
697:         USE KEY, ONLY: R2INTS, DIPINTS 
698:  
699:         IF (ALLOCATED(R2INTS)) DEALLOCATE(R2INTS) 
700:         IF (ALLOCATED(DIPINTS)) DEALLOCATE(DIPINTS) 
701:  
702:     END SUBROUTINE ORBITALS_FINISH 
703:  
704:     SUBROUTINE CHECK_FILE_EXISTS(FNAME, EX) 
705:  
706:         ! Check if a file of the given name exists, exit if it does not and EX=.TRUE. 
707:  
708:         LOGICAL :: FILE_EXISTS, TEX 
709:         LOGICAL, INTENT(IN), OPTIONAL :: EX 
710:         CHARACTER(*), INTENT(IN) :: FNAME 
711:  
712:         IF(PRESENT(EX)) THEN 
713:           TEX = EX 
714:         ELSE 
715:           TEX = .FALSE. 
716:         ENDIF 
717:  
718:         INQUIRE(FILE=FNAME, EXIST=FILE_EXISTS) 
719:         IF(.NOT. FILE_EXISTS) THEN 
720:           WRITE(*,'(A)') "The file " // TRIM(FNAME) // " does not exist." 
721:           IF(TEX.EQV..TRUE.) THEN 
722:             STOP 
723:           ENDIF 
724:         ENDIF 
725:  
726:     END SUBROUTINE CHECK_FILE_EXISTS 
727:  
728: END MODULE ORBITALS_MOD 


r32759/potential.f 2017-06-09 15:30:15.047955081 +0100 r32758/potential.f 2017-06-09 15:30:16.243971716 +0100
 31:          USE FINITE_DIFFERENCES 31:          USE FINITE_DIFFERENCES
 32:          USE MODAMBER9,only : ifswitch,goodstructure1,irespa,cisarray1,checkcistransalways,checkcistransalwaysdna, 32:          USE MODAMBER9,only : ifswitch,goodstructure1,irespa,cisarray1,checkcistransalways,checkcistransalwaysdna,
 33:      1   checkcistransalwaysrna 33:      1   checkcistransalwaysrna
 34:          ! hk286 34:          ! hk286
 35:          USE GENRIGID 35:          USE GENRIGID
 36:          USE MULTIPOT, ONLY: MULTIPOT_CALL 36:          USE MULTIPOT, ONLY: MULTIPOT_CALL
 37:          USE CHIRALITY, ONLY: CIS_TRANS_CHECK, CHIRALITY_CHECK 37:          USE CHIRALITY, ONLY: CIS_TRANS_CHECK, CHIRALITY_CHECK
 38:          USE AMBER12_INTERFACE_MOD, ONLY: AMBER12_ENERGY_AND_GRADIENT, POT_ENE_REC_C, AMBER12_NUM_HESS 38:          USE AMBER12_INTERFACE_MOD, ONLY: AMBER12_ENERGY_AND_GRADIENT, POT_ENE_REC_C, AMBER12_NUM_HESS
 39:          USE MODCUDALBFGS, ONLY: CUDA_ENEGRAD_WRAPPER 39:          USE MODCUDALBFGS, ONLY: CUDA_ENEGRAD_WRAPPER
 40:          USE OPEP_INTERFACE_MOD, ONLY: OPEP_ENERGY_AND_GRADIENT,OPEP_NUM_HESS 40:          USE OPEP_INTERFACE_MOD, ONLY: OPEP_ENERGY_AND_GRADIENT,OPEP_NUM_HESS
 41:          USE ORBITALS_MOD, ONLY: GET_ORBITAL_LOCALITY 
 42:          ! use AMHGLOBALS 41:          ! use AMHGLOBALS
 43:  42: 
 44:          IMPLICIT NONE 43:          IMPLICIT NONE
 45:  44: 
 46:          DOUBLE PRECISION COORDS(NOPT) 45:          DOUBLE PRECISION COORDS(NOPT)
 47:          DOUBLE PRECISION ENERGY 46:          DOUBLE PRECISION ENERGY
 48:          DOUBLE PRECISION VNEW(NOPT) 47:          DOUBLE PRECISION VNEW(NOPT)
 49:          LOGICAL GTEST, STEST 48:          LOGICAL GTEST, STEST
 50:          DOUBLE PRECISION RMS 49:          DOUBLE PRECISION RMS
 51:          LOGICAL PTEST, BOXTEST, file_exists 50:          LOGICAL PTEST, BOXTEST, file_exists
549: !                     ENDIF548: !                     ENDIF
550: !                  ENDDO549: !                  ENDDO
551: !                ENDDO550: !                ENDDO
552: !                STOP551: !                STOP
553:             ELSEIF (MLPB3T) THEN552:             ELSEIF (MLPB3T) THEN
554:                CALL MLPB3(COORDS, VNEW, ENERGY, GTEST, SSTEST)553:                CALL MLPB3(COORDS, VNEW, ENERGY, GTEST, SSTEST)
555:             ELSEIF (MLP3T) THEN554:             ELSEIF (MLP3T) THEN
556:                CALL MLP3(COORDS, VNEW, ENERGY, GTEST, SSTEST)555:                CALL MLP3(COORDS, VNEW, ENERGY, GTEST, SSTEST)
557:             ELSEIF (MLQT) THEN556:             ELSEIF (MLQT) THEN
558:                CALL MLQ(COORDS, VNEW, ENERGY, GTEST, SSTEST)557:                CALL MLQ(COORDS, VNEW, ENERGY, GTEST, SSTEST)
559:             ELSEIF (ORBITALS) THEN 
560:                IF (ALLOCATED(HESS)) DEALLOCATE(HESS) 
561:                IF (.NOT.(ALLOCATED(HESS))) THEN 
562:                     PRINT*,'allocating hessian' 
563:                     ALLOCATE(HESS(1:NROTS,1:NROTS)) 
564:                END IF 
565:                CALL GET_ORBITAL_LOCALITY(COORDS, VNEW, ENERGY, GTEST, SSTEST) 
566:                DIFF=1.0D-6 
567:                PRINT*,'analytic and numerical gradients:' 
568:                IF (.NOT.(ALLOCATED(HESS))) THEN 
569:                     PRINT*,'allocating hessian' 
570:                     ALLOCATE(HESS(NROTS,NROTS)) 
571:                END IF 
572:                !CALL MLPVB3(COORDS, VNEW, ENERGY, .TRUE., .TRUE.) 
573:                CALL GET_ORBITAL_LOCALITY(COORDS, VNEW, ENERGY, .TRUE., .TRUE.) 
574:                !PRINT '(A,I8)','SIZE(HESS)=',SIZE(HESS) 
575:                !HESSDUM(1:NROTS,1:NROTS)=HESS(1:NROTS,1:NROTS) 
576:                !DO J1=1,NATOMS 
577:                !   COORDS(J1)=COORDS(J1)+DIFF 
578: !                  CALL MLPVB3(COORDS,VPLUS,EPLUS,.FALSE.,.FALSE.) 
579:                !   CALL GET_ORBITAL_LOCALITY(COORDS, VPLUS, EPLUS, .FALSE., .FALSE.) 
580:                !   COORDS(J1)=COORDS(J1)-2.0D0*DIFF 
581: !                  CALL MLPVB3(COORDS,VMINUS,EMINUS,.FALSE.,.FALSE.) 
582:                !   CALL GET_ORBITAL_LOCALITY(COORDS, VMINUS, EMINUS, .FALSE., .FALSE.) 
583:                !   COORDS(J1)=COORDS(J1)+DIFF 
584:                !   IF ((ABS(VNEW(J1)).NE.0.0D0).AND.  
585:      &         !      (ABS(100.0D0*(VNEW(J1)-(EPLUS-EMINUS)/(2.0D0*DIFF))/VNEW(J1)).GT.1.0D0)) THEN 
586:                !      WRITE(*,'(A,I5,2F20.10,A)') 'anal and num ',J1,VNEW(J1),(EPLUS-EMINUS)/(2.0D0*DIFF),'   X' 
587:                !   ELSE 
588:                !      WRITE(*,'(A,I5,2F20.10)') 'anal and num ',J1,VNEW(J1),(EPLUS-EMINUS)/(2.0D0*DIFF) 
589:                !   ENDIF 
590:                !ENDDO 
591:                !PRINT*,'analytic and numerical second derivatives:' 
592:                !DO J1=1,NATOMS 
593:                   !COORDS(J1)=COORDS(J1)+DIFF 
594: !                  CALL MLPVB3(COORDS,VPLUS,EPLUS,.TRUE.,.FALSE.) 
595:                   !CALL GET_ORBITAL_LOCALITY(COORDS, VPLUS, ENERGY, .TRUE., .FALSE.) 
596:                   !COORDS(J1)=COORDS(J1)-2.0D0*DIFF 
597:                   !CALL MLPVB3(COORDS,VMINUS,EMINUS,.TRUE.,.FALSE.) 
598:                   !CALL GET_ORBITAL_LOCALITY(COORDS, VMINUS, ENERGY, .TRUE., .FALSE.) 
599:                   !COORDS(J1)=COORDS(J1)+DIFF 
600:                   !DO J2=1,NATOMS 
601:                   !   IF ((ABS(HESS(J1,J2)).NE.0.0D0).AND.  
602:      &            !      (ABS(100.0D0*(HESS(J1,J2)-(VPLUS(J2)-VMINUS(J2))/(2.0D0*DIFF))/HESS(J1,J2)).GT.1.0D0)) THEN 
603:                   !   WRITE(*,'(2I5,2G20.10,A)') J1,J2,HESS(J1,J2),(VPLUS(J2)-VMINUS(J2))/(2.0D0*DIFF),'   X' 
604:                   !   ELSE 
605:                   !      WRITE(*,'(2I5,2G20.10,A)') J1,J2,HESS(J1,J2),(VPLUS(J2)-VMINUS(J2))/(2.0D0*DIFF) 
606:                   !   ENDIF 
607:                   !ENDDO 
608:                 !ENDDO 
609:                !STOP 
610:             ELSE558:             ELSE
611:                CALL FUNCTIONAL( COORDS, VNEW, ENERGY, GTEST, SSTEST)559:                CALL FUNCTIONAL( COORDS, VNEW, ENERGY, GTEST, SSTEST)
612:             ENDIF560:             ENDIF
613:             IF (PTEST.OR.MLPPROB.OR.MLQPROB) THEN561:             IF (PTEST.OR.MLPPROB.OR.MLQPROB) THEN
614:                 WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,'         '562:                 WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,'         '
615:                 WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,'         '563:                 WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,'         '
616:             ENDIF564:             ENDIF
617: !           IF (MLPPROB) STOP565: !           IF (MLPPROB) STOP
618:             ! CALL CTEST(NATOMS, COORDS, VNEW, ENERGY, GTEST, SSTEST)566:             ! CALL CTEST(NATOMS, COORDS, VNEW, ENERGY, GTEST, SSTEST)
619:             ! CALL TWODFUNC(COORDS,VNEW,ENERGY,GTEST,SSTEST)567:             ! CALL TWODFUNC(COORDS,VNEW,ENERGY,GTEST,SSTEST)


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0