hdiff output

r30266/key.f90 2016-04-14 11:30:13.529653526 +0100 r30265/key.f90 2016-04-14 11:30:14.297664068 +0100
 45:      &        QSPCFWT, QTIP4PFT, CFUSIONT, DUMPINTXYZ, DUMPINTEOS, INTLJT, INTTST, EYTRAPT, OHCELLT, MKTRAPT, & 45:      &        QSPCFWT, QTIP4PFT, CFUSIONT, DUMPINTXYZ, DUMPINTEOS, INTLJT, INTTST, EYTRAPT, OHCELLT, MKTRAPT, &
 46:      &        INTFREEZET, LPERMDIST, CHECKNEGATIVET, CHECKOVERLAPT, ACK1, ACK2, CONDATT, USERPOTT, & 46:      &        INTFREEZET, LPERMDIST, CHECKNEGATIVET, CHECKOVERLAPT, ACK1, ACK2, CONDATT, USERPOTT, &
 47:      &        CONCUTFRACT, CONCUTABST, ENDNUMHESS2, CHARMMDFTBT, PAIRCOLOURT, REVERSEUPHILLT, WHOLEDNEB, & 47:      &        CONCUTFRACT, CONCUTABST, ENDNUMHESS2, CHARMMDFTBT, PAIRCOLOURT, REVERSEUPHILLT, WHOLEDNEB, &
 48:      &        NONEBMAX, READMASST, ONEDAPBCT, ONEDPBCT, INVTONEDPBCT, INVTTWODPBCT, TWODAPBCT, TWODPBCT, THREEDAPBCT, & 48:      &        NONEBMAX, READMASST, ONEDAPBCT, ONEDPBCT, INVTONEDPBCT, INVTTWODPBCT, TWODAPBCT, TWODPBCT, THREEDAPBCT, &
 49:      &        THREEDPBCT, FOURDAPBCT, FOURDPBCT, MODEDOWNT, CHEMSHIFT, TTM3T, & 49:      &        THREEDPBCT, FOURDAPBCT, FOURDPBCT, MODEDOWNT, CHEMSHIFT, TTM3T, &
 50:      &        NOINVERSION, INVERTPT, KNOWVECS, PMPATHT, AAORIENTT, MULTIJOBT, QUIPARGSTRT, QUIPPARAMST, HESSDUMPT, & 50:      &        NOINVERSION, INVERTPT, KNOWVECS, PMPATHT, AAORIENTT, MULTIJOBT, QUIPARGSTRT, QUIPPARAMST, HESSDUMPT, &
 51:      &        CLASSICALRATEST, TSPLITTINGT, HESSREADT, INSTANTONOPTT,INSTANTONSTARTDUMPT,VARSTEPOPTT, MOLPRO, REAXFFT, & 51:      &        CLASSICALRATEST, TSPLITTINGT, HESSREADT, INSTANTONOPTT,INSTANTONSTARTDUMPT,VARSTEPOPTT, MOLPRO, REAXFFT, &
 52:      &        EIGENONLY,OVERCONV, GLJT,CLSTRINGT,CLSTRINGTST, PHI4MODT, EX1DT, MCPATHT, MCBIAST, RPHT, TWISTT, MCPATH2T, & 52:      &        EIGENONLY,OVERCONV, GLJT,CLSTRINGT,CLSTRINGTST, PHI4MODT, EX1DT, MCPATHT, MCBIAST, RPHT, TWISTT, MCPATH2T, &
 53:      &        PBST, SSHT, GAUSSIAN03, CPPNEBT, CUDAT, CUDATIMET, TRUSTMODET,MODELOST, METRICTENSOR, INTSPRINGACTIVET, & 53:      &        PBST, SSHT, GAUSSIAN03, CPPNEBT, CUDAT, CUDATIMET, TRUSTMODET,MODELOST, METRICTENSOR, INTSPRINGACTIVET, &
 54:      &        PERMGUESS, QCIPERMCHECK, DUMPFRQST, MULTIPOTT, MLP3T, MLPB3T, DUMPBESTPATH, ALIGNRBST, AVOID_COLLISIONS, MLPPROB, & 54:      &        PERMGUESS, QCIPERMCHECK, DUMPFRQST, MULTIPOTT, MLP3T, MLPB3T, DUMPBESTPATH, ALIGNRBST, AVOID_COLLISIONS, MLPPROB, &
 55:      &        MALONALDEHYDE, MLPNEWREG, DJWRBT, STEALTHYT 55:      &        MALONALDEHYDE, MLPNEWREG, DJWRBT
 56:  56: 
 57: ! bf269 > polymer in a pore (non-bonding (LJ) energy from neighbours is not subtracted) 57: ! bf269 > polymer in a pore (non-bonding (LJ) energy from neighbours is not subtracted)
 58:       LOGICAL :: PORE8T = .FALSE. ! add 8th power cylindrical pore to the potential? 58:       LOGICAL :: PORE8T = .FALSE. ! add 8th power cylindrical pore to the potential?
 59:       INTEGER :: PORE8_AXIS = 3 ! principal axis of the cylindric pore (1:x, 2:y, 3:z) 59:       INTEGER :: PORE8_AXIS = 3 ! principal axis of the cylindric pore (1:x, 2:y, 3:z)
 60:       DOUBLE PRECISION :: PORE8_ENERGY = 1.0d1 ! energy of the pore when radius = 1 60:       DOUBLE PRECISION :: PORE8_ENERGY = 1.0d1 ! energy of the pore when radius = 1
 61:       LOGICAL :: HARMPOLYT = .FALSE. ! add harmonic bonds between the beads 61:       LOGICAL :: HARMPOLYT = .FALSE. ! add harmonic bonds between the beads
 62:       DOUBLE PRECISION :: HARMPOLY_BONLEN = 0.0d0 ! equilibrium length of springs between beads 62:       DOUBLE PRECISION :: HARMPOLY_BONLEN = 0.0d0 ! equilibrium length of springs between beads
 63:       DOUBLE PRECISION :: HARMPOLY_K = 1.0d2 ! force constant of the springs 63:       DOUBLE PRECISION :: HARMPOLY_K = 1.0d2 ! force constant of the springs
 64:  64: 
 65: ! hk286 > generalised THOMSON problem 65: ! hk286 > generalised THOMSON problem
 91:      &        BISECTMAXENERGY, BISECTMINDIST, BLFACTOR, NEBRESEEDEMAX, NEBRESEEDBMAX, NEBRESEEDDEL1, & 91:      &        BISECTMAXENERGY, BISECTMINDIST, BLFACTOR, NEBRESEEDEMAX, NEBRESEEDBMAX, NEBRESEEDDEL1, &
 92:      &        NEBRESEEDDEL2, INTCONSTRAINTTOL, INTCONSTRAINTDEL, RBCUTOFF, INTCONSTRAINTREP, INTCONSTRAINREPCUT, & 92:      &        NEBRESEEDDEL2, INTCONSTRAINTTOL, INTCONSTRAINTDEL, RBCUTOFF, INTCONSTRAINTREP, INTCONSTRAINREPCUT, &
 93:      &        REDOK, REDOFRAC, D1INIT, D2INIT, REDOE1, REDOE2, RPBETA, REPCON, PFORCE, & 93:      &        REDOK, REDOFRAC, D1INIT, D2INIT, REDOE1, REDOE2, RPBETA, REPCON, PFORCE, &
 94:      &        CPCONSTRAINTTOL, CPCONSTRAINTDEL, CPCONSTRAINTREP, CPCONSTRAINREPCUT, CPCONFRAC, & 94:      &        CPCONSTRAINTTOL, CPCONSTRAINTDEL, CPCONSTRAINTREP, CPCONSTRAINREPCUT, CPCONFRAC, &
 95:      &        INTLJTOL, INTLJDEL, INTLJEPS, IMSEPMIN, IMSEPMAX, TRAPK, MINOVERLAP, & 95:      &        INTLJTOL, INTLJDEL, INTLJEPS, IMSEPMIN, IMSEPMAX, TRAPK, MINOVERLAP, &
 96:      &        INTFREEZETOL, LOCALPERMCUT, LOCALPERMCUT2, LOCALPERMCUTINC, CHECKREPCUTOFF, CONCUTABS, & 96:      &        INTFREEZETOL, LOCALPERMCUT, LOCALPERMCUT2, LOCALPERMCUTINC, CHECKREPCUTOFF, CONCUTABS, &
 97:      &        CONCUTFRAC, ENDNUMHESSDELTA, DNEBEFRAC, QCHEMSCALE, KAA, SIGMAAA, QUIPATOMMASS, TEMPERATURE1, & 97:      &        CONCUTFRAC, ENDNUMHESSDELTA, DNEBEFRAC, QCHEMSCALE, KAA, SIGMAAA, QUIPATOMMASS, TEMPERATURE1, &
 98:      &        DISTORTINST,DELTAINST,MOLPROSCALE,COVER,STTSRMSCONV,LAN_DIST,LANCONV,LANFACTOR, & 98:      &        DISTORTINST,DELTAINST,MOLPROSCALE,COVER,STTSRMSCONV,LAN_DIST,LANCONV,LANFACTOR, &
 99:      &        STOCKEXP, JPARAM, MCPATHTEMP, MCPATHDMAX, MCPATHSTEP, MCPATHACCRATIO, BIASFAC, & 99:      &        STOCKEXP, JPARAM, MCPATHTEMP, MCPATHDMAX, MCPATHSTEP, MCPATHACCRATIO, BIASFAC, &
100:      &        MCADDDEV, MCPATHQMIN, MCPATHQMAX, RPHQMIN, RPHQMAX, RPHTEMP, TWISTF, TWISTREF, MCPATHADDREF, &100:      &        MCADDDEV, MCPATHQMIN, MCPATHQMAX, RPHQMIN, RPHQMAX, RPHTEMP, TWISTF, TWISTREF, MCPATHADDREF, &
101:      &        MCPATHGWS, MCPATHGWQ, MCPATHNEGLECT, MCPATHTOL, FRAMESDIFF,TMRATIO, INTMINFAC, MLPLAMBDA, COLL_TOL, KLIM, SCA101:      &        MCPATHGWS, MCPATHGWQ, MCPATHNEGLECT, MCPATHTOL, FRAMESDIFF,TMRATIO, INTMINFAC, MLPLAMBDA, COLL_TOL
102: 102: 
103: !     sf344103: !     sf344
104:       DOUBLE PRECISION :: PCUTOFF,PYA11(3),PYA21(3),PYA12(3),PYA22(3),PEPSILON1(3),PSCALEFAC1(2),PSCALEFAC2(2), &104:       DOUBLE PRECISION :: PCUTOFF,PYA11(3),PYA21(3),PYA12(3),PYA22(3),PEPSILON1(3),PSCALEFAC1(2),PSCALEFAC2(2), &
105:      &                     PEPSILONATTR(2),PSIGMAATTR(2), PYOVERLAPTHRESH, LJSITECOORDS(3), LJGSITESIGMA, LJGSITEEPS, &105:      &                     PEPSILONATTR(2),PSIGMAATTR(2), PYOVERLAPTHRESH, LJSITECOORDS(3), LJGSITESIGMA, LJGSITEEPS, &
106:      &                     PYLOCALSTEP(2)106:      &                     PYLOCALSTEP(2)
107:  107:  
108:       DOUBLE PRECISION, ALLOCATABLE :: POINTSDECA(:), POINTSICOS(:)108:       DOUBLE PRECISION, ALLOCATABLE :: POINTSDECA(:), POINTSICOS(:)
109:       DOUBLE PRECISION, ALLOCATABLE :: VT(:), pya1bin(:,:),pya2bin(:,:)109:       DOUBLE PRECISION, ALLOCATABLE :: VT(:), pya1bin(:,:),pya2bin(:,:)
110:       LOGICAL          :: LJSITE,BLJSITE,LJSITEATTR,PYBINARYT,PARAMONOVPBCX,PARAMONOVPBCY,PARAMONOVPBCZ,PARAMONOVCUTOFF110:       LOGICAL          :: LJSITE,BLJSITE,LJSITEATTR,PYBINARYT,PARAMONOVPBCX,PARAMONOVPBCY,PARAMONOVPBCZ,PARAMONOVCUTOFF
111:       LOGICAL          :: PYGPERIODICT,ELLIPSOIDT,LJSITECOORDST,REALIGNXYZ,MULTISITEPYT,LJGSITET,NORMALMODET111:       LOGICAL          :: PYGPERIODICT,ELLIPSOIDT,LJSITECOORDST,REALIGNXYZ,MULTISITEPYT,LJGSITET,NORMALMODET


r30266/keywords.f 2016-04-14 11:30:13.725656223 +0100 r30265/keywords.f 2016-04-14 11:30:14.485666648 +0100
5615: ! SQVV allows the first NIterSQVVGuessMax DNEB iterations to be done using5615: ! SQVV allows the first NIterSQVVGuessMax DNEB iterations to be done using
5616: ! SQVV - switches to LBFGS minimisation after NIterSQVVGuessMax iterations5616: ! SQVV - switches to LBFGS minimisation after NIterSQVVGuessMax iterations
5617: ! or if the RMS force goes below SQVVGuessRMSTol.5617: ! or if the RMS force goes below SQVVGuessRMSTol.
5618: ! 5618: ! 
5619:          ELSE IF (WORD == 'SQVV') THEN5619:          ELSE IF (WORD == 'SQVV') THEN
5620:             SQVVGUESS=.TRUE.5620:             SQVVGUESS=.TRUE.
5621:             IF (NITEMS.GT.1) CALL READI(NITERSQVVGUESSMAX)5621:             IF (NITEMS.GT.1) CALL READI(NITERSQVVGUESSMAX)
5622:             IF (NITEMS.GT.2) CALL READF(SQVVGUESSRMSTOL)5622:             IF (NITEMS.GT.2) CALL READF(SQVVGUESSRMSTOL)
5623:          ELSE IF (WORD.EQ.'SSH') THEN5623:          ELSE IF (WORD.EQ.'SSH') THEN
5624:             SSHT=.TRUE.5624:             SSHT=.TRUE.
5625:  
5626:          ELSE IF (WORD.EQ.'STEALTHY') THEN 
5627:             STEALTHYT=.TRUE. 
5628:             CALL READF(KLIM) 
5629:             IF (NITEMS.GT.2) THEN 
5630:                CALL READF(SCA) 
5631:             ELSE 
5632:                SCA=1 
5633:             END IF 
5634: ! 5625: ! 
5635: ! NSTEPMIN sets the minimum number of steps allowed before convergence.5626: ! NSTEPMIN sets the minimum number of steps allowed before convergence.
5636: ! 5627: ! 
5637:          ELSE IF (WORD .EQ. 'STEPMIN') THEN5628:          ELSE IF (WORD .EQ. 'STEPMIN') THEN
5638:             CALL READI(NSTEPMIN)5629:             CALL READI(NSTEPMIN)
5639: ! 5630: ! 
5640: ! STEPS n sets the number of optimisation steps to perform5631: ! STEPS n sets the number of optimisation steps to perform
5641: ! per call to OPTIM                                    - default n=15632: ! per call to OPTIM                                    - default n=1
5642: ! If BFGSSTEPS is not specified then it is set to the same value as NSTEPS5633: ! If BFGSSTEPS is not specified then it is set to the same value as NSTEPS
5643: ! 5634: ! 


r30266/potential.f 2016-04-14 11:30:13.917658851 +0100 r30265/potential.f 2016-04-14 11:30:14.681669342 +0100
3499:                WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY3499:                WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY
3500:             ENDIF3500:             ENDIF
3501:             ! davidg: introduced userpot here3501:             ! davidg: introduced userpot here
3502:          ELSE IF (USERPOTT) THEN3502:          ELSE IF (USERPOTT) THEN
3503:             CALL USERPOT_POTENTIALHESS(3*NATOMS,COORDS,VNEW,ENERGY,GTEST,SSTEST,HESS)3503:             CALL USERPOT_POTENTIALHESS(3*NATOMS,COORDS,VNEW,ENERGY,GTEST,SSTEST,HESS)
3504:             IF (PTEST) THEN3504:             IF (PTEST) THEN
3505:                WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,' meV'3505:                WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,' meV'
3506:                WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,' meV'3506:                WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,' meV'
3507:             ENDIF3507:             ENDIF
3508: 3508: 
3509:             ELSE IF (STEALTHYT) THEN 
3510:                CALL STEALTHY(COORDS, VNEW, ENERGY, GTEST, SSTEST) 
3511:   
3512:                !DIFF=1.0D-6 
3513:                !PRINT*,'ANALYTIC AND NUMERICAL GRADIENTS:NATOMS=',NATOMS 
3514:                !DO J1=1, 3*NATOMS 
3515:                !   COORDS(J1)=COORDS(J1)+DIFF 
3516:                !   CALL STEALTHY(COORDS,VPLUS,EPLUS,.FALSE.,.FALSE.) 
3517:                !   COORDS(J1)=COORDS(J1)-2.0D0*DIFF 
3518:                !   CALL STEALTHY(COORDS,VMINUS,EMINUS,.FALSE.,.FALSE.) 
3519:                !   COORDS(J1)=COORDS(J1)+DIFF 
3520:                ! 
3521:                !  IF ((ABS(VNEW(J1)).NE.0.0D0).AND.(100.0D0*ABS((VNEW(J1)-(EPLUS-EMINUS)/(2.0D0*DIFF))/VNEW(J1)).GT.1.0D-3)) THEN 
3522:                !      WRITE(*,'(I5,2G20.10,A)') J1,VNEW(J1),(EPLUS-EMINUS)/(2.0D0*DIFF), '  X' 
3523:                !   ELSE 
3524:                !      WRITE(*,'(I5,2G20.10,A)') J1,VNEW(J1),(EPLUS-EMINUS)/(2.0D0*DIFF) 
3525:                !   END IF 
3526:                !END DO 
3527:    
3528:                !CALL STEALTHY(COORDS, VNEW, ENERGY, .TRUE., .TRUE.) 
3529:                !OPEN(UNIT=1499, FILE="hessian-matrix") 
3530:                !WRITE(1499,*),'ANALYTIC AND NUMERICAL SECOND DERIVATIVES:' 
3531:                !DO J1=1, 3*NATOMS 
3532:                !   COORDS(J1)=COORDS(J1)+DIFF 
3533:                !   CALL STEALTHY(COORDS,VPLUS,EPLUS,.TRUE.,.FALSE.) 
3534:                !   COORDS(J1)=COORDS(J1)-2.0D0*DIFF 
3535:                !   CALL STEALTHY(COORDS,VMINUS,EMINUS,.TRUE.,.FALSE.) 
3536:                !   COORDS(J1)=COORDS(J1)+DIFF 
3537:                !   DO J2=1, 3*NATOMS 
3538:                !      IF (ABS(HESS(J1,J2)-(VPLUS(J2)-VMINUS(J2))/(2.0D0*DIFF)).GT.(1.0D-5*HESS(J1,J2))) THEN 
3539:                !         WRITE(1499,'(2I5,2F20.15,A)') J1,J2,HESS(J1,J2),(VPLUS(J2)-VMINUS(J2))/(2.0D0*DIFF),'  X' 
3540:                !      ELSE 
3541:                !         WRITE(1499,'(2I5,2F20.15)') J1,J2,HESS(J1,J2),(VPLUS(J2)-VMINUS(J2))/(2.0D0*DIFF) 
3542:                !      END IF 
3543:                !   END DO 
3544:                !END DO 
3545:  
3546:  
3547: 3509: 
3548: 3510: 
3549:             ! 3511:             ! 
3550:             ! Check numerical first and second derivatives3512:             ! Check numerical first and second derivatives
3551:             ! 3513:             ! 
3552:             ! DIFF=1.0D-43514:             ! DIFF=1.0D-4
3553:             ! PRINT*,'analytic and numerical gradients:'3515:             ! PRINT*,'analytic and numerical gradients:'
3554:             ! DO J1=1,3*NATOMS3516:             ! DO J1=1,3*NATOMS
3555:             ! COORDS(J1)=COORDS(J1)+DIFF3517:             ! COORDS(J1)=COORDS(J1)+DIFF
3556:             ! CALL FCAPSID(NATOMS,COORDS,HDUMM,VPLUS,EPLUS,.FALSE.,.FALSE.)3518:             ! CALL FCAPSID(NATOMS,COORDS,HDUMM,VPLUS,EPLUS,.FALSE.,.FALSE.)


r30266/stealthy.f90 2016-04-14 11:30:14.101661377 +0100 r30265/stealthy.f90 2016-04-14 11:30:14.873671977 +0100
  1: SUBROUTINE STEALTHY(XS,TIDU,ENERGY,GTEST,STEST)  1: svn: E195012: Unable to find repository location for 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/OPTIM/source/stealthy.f90' in revision 30265
  2:     USE COMMONS 
  3:     USE KEY 
  4:     USE MODHESS 
  5:     IMPLICIT NONE 
  6:     LOGICAL GTEST,STEST 
  7:     INTEGER M, M1, M2, M3, J1, J2, I1, I2, L1, L2, KI 
  8:     DOUBLE PRECISION XS(3*NATOMS), KX(3), KY(3), KZ(3), KN(3), TIDU(3*NATOMS), & 
  9:                      CKR(NATOMS), SKR(NATOMS) 
 10:     DOUBLE PRECISION PI, KR ,C, S, ENERGY, VF, IM, KM 
 11:  
 12:     IF (.NOT.ALLOCATED(HESS)) ALLOCATE(HESS(3*NATOMS, 3*NATOMS)) 
 13:  
 14: !   initialize 
 15:     PI=3.141592653589793 
 16:     KM=KLIM*KLIM 
 17:     KI=DINT(KLIM) 
 18:     ENERGY=0 
 19:  
 20:     IF (GTEST) THEN 
 21:        TIDU=0 
 22:     END IF 
 23:  
 24:     IF (STEST) THEN 
 25:        HESS=0 
 26:     END IF 
 27:  
 28: !   get the initial box 
 29:  
 30: !   calculate the wave vector 
 31:     KX(1)=2*PI/BULK_BOXVEC(1) 
 32:     KX(2)=0 
 33:     KX(3)=0 
 34:     KY(1)=0 
 35:     KY(2)=2*PI/BULK_BOXVEC(2) 
 36:     KY(3)=0 
 37:     KZ(1)=0 
 38:     KZ(2)=0 
 39:     KZ(3)=2*PI/BULK_BOXVEC(3) 
 40:     VF=BULK_BOXVEC(1)*BULK_BOXVEC(2)*BULK_BOXVEC(3) 
 41: !   WRITE(*,*) KI(1), KI(2), KI(3) 
 42:  
 43: !   periodic boundary condition 
 44:     IF (BULKT) THEN 
 45:        DO J1=1, NATOMS 
 46:           J2=3*(J1-1) 
 47:           XS(J2+1)=XS(J2+1) - BULK_BOXVEC(1)*DNINT(XS(J2+1)/BULK_BOXVEC(1)) 
 48:           XS(J2+2)=XS(J2+2) - BULK_BOXVEC(2)*DNINT(XS(J2+2)/BULK_BOXVEC(2)) 
 49:           XS(J2+3)=XS(J2+3) - BULK_BOXVEC(3)*DNINT(XS(J2+3)/BULK_BOXVEC(3)) 
 50:        END DO 
 51:     END IF 
 52:   
 53: !   sum of k 
 54:     DO M1=0, KI 
 55:        DO M2=-KI, KI 
 56:           DO M3=-KI, KI 
 57:  
 58:              IF (M1==0.AND.M2>=0.AND.M3<0) CYCLE 
 59:              IF (M1==0.AND.M2<=0.AND.M3<=0) CYCLE 
 60:  
 61:              M=M1*M1+M2*M2+M3*M3 
 62:              IF (M>KM.OR.M==0) CYCLE 
 63:  
 64:              C=0 
 65:              S=0 
 66:              KN(1)=KX(1)*M1 + KY(1)*M2 + KZ(1)*M3 
 67:              KN(2)=KX(2)*M1 + KY(2)*M2 + KZ(2)*M3 
 68:              KN(3)=KX(3)*M1 + KY(3)*M2 + KZ(3)*M3 
 69:  
 70: !            calculation of Re[n(k)] and Im[n(k)] 
 71:              DO I1=1, NATOMS 
 72:                 I2=3*(I1-1) 
 73:                 KR=KN(1)*XS(I2+1) + KN(2)*XS(I2+2) + KN(3)*XS(I2+3) 
 74:                 CKR(I1)=COS(KR) 
 75:                 SKR(I1)=SIN(KR) 
 76:                 C=C+CKR(I1) 
 77:                 S=S+SKR(I1) 
 78:              END DO 
 79:  
 80: !            E=n(k)^2/Vf 
 81:              ENERGY=ENERGY+ SCA*(C*C+S*S)/VF 
 82:  
 83: !            gradient 
 84:              IF (GTEST) THEN 
 85:                 DO L1=1, NATOMS 
 86:                    L2=3*(L1-1) 
 87:                    IM=( C*SKR(L1) - S*CKR(L1) )*SCA 
 88:                    TIDU(L2+1)=TIDU(L2+1) - 2*KN(1)*IM/VF 
 89:                    TIDU(L2+2)=TIDU(L2+2) - 2*KN(2)*IM/VF 
 90:                    TIDU(L2+3)=TIDU(L2+3) - 2*KN(3)*IM/VF 
 91:                 END DO 
 92:              END IF 
 93:  
 94: !            hessian matrix 
 95:              IF (STEST) THEN 
 96:                 CALL STHESS(XS, KN, VF, CKR, C, SKR, S) 
 97:              END IF 
 98:  
 99:           END DO 
100:        END DO 
101:     END DO 
102:     
103: !   change the box back to cubic 
104:  
105: END SUBROUTINE 
106:  
107: SUBROUTINE STHESS(X, K, V, CK, CSUM, SK, SSUM) 
108:    USE COMMONS 
109:    USE KEY 
110:    USE MODHESS 
111:    IMPLICIT NONE 
112:    INTEGER I1, I2, J1, J2 
113:    DOUBLE PRECISION X(3*NATOMS), K(3), RIJ(3), CK(NATOMS), SK(NATOMS) 
114:    DOUBLE PRECISION KRIJ, V, CSUM, SSUM, CI, SI, CF, CSKR 
115:  
116:    DO I1=1, NATOMS 
117:       I2=3*(I1-1) 
118:  
119:       CI=CSUM-CK(I1) 
120:       SI=SSUM-SK(I1) 
121:       CF=(CK(I1)*CI + SK(I1)*SI)*SCA 
122:  
123:       HESS(I2+1, I2+1)=HESS(I2+1, I2+1) - 2*K(1)*K(1)/V * CF 
124:       HESS(I2+1, I2+2)=HESS(I2+1, I2+2) - 2*K(1)*K(2)/V * CF 
125:       HESS(I2+1, I2+3)=HESS(I2+1, I2+3) - 2*K(1)*K(3)/V * CF 
126:       HESS(I2+2, I2+2)=HESS(I2+2, I2+2) - 2*K(2)*K(2)/V * CF 
127:       HESS(I2+2, I2+3)=HESS(I2+2, I2+3) - 2*K(2)*K(3)/V * CF 
128:       HESS(I2+3, I2+3)=HESS(I2+3, I2+3) - 2*K(3)*K(3)/V * CF 
129:  
130:       DO J1=I1+1, NATOMS 
131:          J2=3*(J1-1) 
132:  
133:          RIJ(1)=X(I2+1) - X(J2+1) 
134:          RIJ(2)=X(I2+2) - X(J2+2) 
135:          RIJ(3)=X(I2+3) - X(J2+3) 
136:  
137:          KRIJ=K(1)*RIJ(1) + K(2)*RIJ(2) + K(3)*RIJ(3) 
138:          CSKR= COS(KRIJ)*SCA 
139:  
140:          HESS(I2+1, J2+1)=HESS(I2+1, J2+1) + 2*K(1)*K(1)/V * CSKR 
141:          HESS(I2+1, J2+2)=HESS(I2+1, J2+2) + 2*K(1)*K(2)/V * CSKR 
142:          HESS(I2+1, J2+3)=HESS(I2+1, J2+3) + 2*K(1)*K(3)/V * CSKR 
143:  
144:          HESS(I2+2, J2+1)=HESS(I2+2, J2+1) + 2*K(2)*K(1)/V * CSKR 
145:          HESS(I2+2, J2+2)=HESS(I2+2, J2+2) + 2*K(2)*K(2)/V * CSKR 
146:          HESS(I2+2, J2+3)=HESS(I2+2, J2+3) + 2*K(2)*K(3)/V * CSKR 
147:  
148:          HESS(I2+3, J2+1)=HESS(I2+3, J2+1) + 2*K(3)*K(1)/V * CSKR 
149:          HESS(I2+3, J2+2)=HESS(I2+3, J2+2) + 2*K(3)*K(2)/V * CSKR 
150:          HESS(I2+3, J2+3)=HESS(I2+3, J2+3) + 2*K(3)*K(3)/V * CSKR 
151:  
152:       END DO 
153:    END DO 
154:  
155:    DO I1=1, NATOMS 
156:       I2=3*(I1-1) 
157:  
158:       HESS(I2+2, I2+1)=HESS(I2+1, I2+2) 
159:       HESS(I2+3, I2+1)=HESS(I2+1, I2+3) 
160:       HESS(I2+3, I2+2)=HESS(I2+2, I2+3) 
161:  
162:       DO J1=I1+1, NATOMS 
163:          J2=3*(J1-1) 
164:  
165:          HESS(J2+1, I2+1)=HESS(I2+1, J2+1) 
166:          HESS(J2+2, I2+1)=HESS(I2+1, J2+2) 
167:          HESS(J2+3, I2+1)=HESS(I2+1, J2+3) 
168:  
169:          HESS(J2+1, I2+2)=HESS(I2+2, J2+1) 
170:          HESS(J2+2, I2+2)=HESS(I2+2, J2+2) 
171:          HESS(J2+3, I2+2)=HESS(I2+2, J2+3) 
172:  
173:          HESS(J2+1, I2+3)=HESS(I2+3, J2+1) 
174:          HESS(J2+2, I2+3)=HESS(I2+3, J2+2) 
175:          HESS(J2+3, I2+3)=HESS(I2+3, J2+3) 
176:  
177:       END DO 
178:    END DO 
179:  
180: END SUBROUTINE 


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0