hdiff output

r31876/commons.f90 2017-02-03 06:30:31.415234392 +0000 r31875/commons.f90 2017-02-03 06:30:32.663251185 +0000
112:      &        CHARMMDFTBT, PERMINVOPT, BLOCKMOVET, MAXERISE_SET, PYT, BINARY_EXAB, CHIROT, POLYT, SANDBOXT, &112:      &        CHARMMDFTBT, PERMINVOPT, BLOCKMOVET, MAXERISE_SET, PYT, BINARY_EXAB, CHIROT, POLYT, SANDBOXT, &
113:      &        RESERVOIRT, DISTOPT, ONEDAPBCT, ONEDPBCT, TWODAPBCT, TWODPBCT, THREEDAPBCT, THREEDPBCT, RATIOT, &113:      &        RESERVOIRT, DISTOPT, ONEDAPBCT, ONEDPBCT, TWODAPBCT, TWODPBCT, THREEDAPBCT, THREEDPBCT, RATIOT, &
114:      &        PTRANDOM, PTINTERVAL, PTSINGLE, PTSETS, CHEMSHIFT, CHEMSHIFT2, CSH, DEBUGss2029, UNIFORMMOVE, RANSEEDT, &114:      &        PTRANDOM, PTINTERVAL, PTSINGLE, PTSETS, CHEMSHIFT, CHEMSHIFT2, CSH, DEBUGss2029, UNIFORMMOVE, RANSEEDT, &
115:      &        TTM3T, NOINVERSION, RIGIDCONTOURT, UPDATERIGIDREFT, HYBRIDMINT, COMPRESSRIGIDT, MWFILMT, &115:      &        TTM3T, NOINVERSION, RIGIDCONTOURT, UPDATERIGIDREFT, HYBRIDMINT, COMPRESSRIGIDT, MWFILMT, &
116:      &        SUPPRESST, MFETT, POLIRT, QUIPT, SWPOTT, MWPOTT, REPMATCHT, GLJT, MLJT, READMASST, SPECMASST, NEWTSALLIST, &116:      &        SUPPRESST, MFETT, POLIRT, QUIPT, SWPOTT, MWPOTT, REPMATCHT, GLJT, MLJT, READMASST, SPECMASST, NEWTSALLIST, &
117:      &        PHI4MODELT, CUDAT, CUDATIMET, AMBER12T, ENERGY_DECOMPT, NEWMOVEST, DUMPMINT, MBPOLT, MOLECULART, GCBHT, SEMIGRAND_MUT, USEROT, & 117:      &        PHI4MODELT, CUDAT, CUDATIMET, AMBER12T, ENERGY_DECOMPT, NEWMOVEST, DUMPMINT, MBPOLT, MOLECULART, GCBHT, SEMIGRAND_MUT, USEROT, & 
118:      &        SAVEMULTIMINONLY, GRADPROBLEMT, INTLJT, CONDATT, QCIPERMCHECK, &118:      &        SAVEMULTIMINONLY, GRADPROBLEMT, INTLJT, CONDATT, QCIPERMCHECK, &
119:      &        INTCONSTRAINTT, INTFREEZET, CHECKCONINT, CONCUTABST, CONCUTFRACT, INTERPCOSTFUNCTION, &119:      &        INTCONSTRAINTT, INTFREEZET, CHECKCONINT, CONCUTABST, CONCUTFRACT, INTERPCOSTFUNCTION, &
120:      &        RBAAT, FREEZENODEST, DUMPINTEOS, DUMPINTXYZ, QCIPOTT, QCIPOT2T, INTSPRINGACTIVET, LPERMDIST, LOCALPERMDIST, QCIRADSHIFTT, &120:      &        RBAAT, FREEZENODEST, DUMPINTEOS, DUMPINTXYZ, QCIPOTT, QCIPOT2T, INTSPRINGACTIVET, LPERMDIST, LOCALPERMDIST, QCIRADSHIFTT, &
121:      &        MLP3T, MKTRAPT, MLPB3T, MLPB3NEWT, MULTIPOTT, QCIAMBERT, MLPNEWREG, DJWRBT, STEALTHYT, LJADDT, QCINOREPINT, RIGIDMDT, &121:      &        MLP3T, MKTRAPT, MLPB3T, MLPB3NEWT, MULTIPOTT, QCIAMBERT, MLPNEWREG, DJWRBT, STEALTHYT, LJADDT, QCINOREPINT, RIGIDMDT, &
122:      &        DUMPMQT, MLQT, MLQPROB, LJADD2T, MLPVB3T, NOREGBIAS, PYADDT122:      &        DUMPMQT, MLQT, MLQPROB, LJADD2T, MLPVB3T, NOREGBIAS
123: !123: !
124:       DOUBLE PRECISION, ALLOCATABLE :: SEMIGRAND_MU(:) 124:       DOUBLE PRECISION, ALLOCATABLE :: SEMIGRAND_MU(:) 
125:       DOUBLE PRECISION, ALLOCATABLE :: ATMASS(:)125:       DOUBLE PRECISION, ALLOCATABLE :: ATMASS(:)
126:       DOUBLE PRECISION, ALLOCATABLE :: SPECMASS(:) 126:       DOUBLE PRECISION, ALLOCATABLE :: SPECMASS(:) 
127: 127: 
128: ! csw34> FREEZEGROUP variables128: ! csw34> FREEZEGROUP variables
129: !129: !
130:       INTEGER :: GROUPCENTRE130:       INTEGER :: GROUPCENTRE
131:       DOUBLE PRECISION :: GROUPRADIUS131:       DOUBLE PRECISION :: GROUPRADIUS
132:       CHARACTER (LEN=2) :: FREEZEGROUPTYPE132:       CHARACTER (LEN=2) :: FREEZEGROUPTYPE
619:       INTEGER :: INTCONMAX=10, NCPCONMAX=10619:       INTEGER :: INTCONMAX=10, NCPCONMAX=10
620:       DOUBLE PRECISION, ALLOCATABLE :: EPSALPHA(:), DISTREF(:)620:       DOUBLE PRECISION, ALLOCATABLE :: EPSALPHA(:), DISTREF(:)
621:       DOUBLE PRECISION, ALLOCATABLE :: CPDISTREF(:)621:       DOUBLE PRECISION, ALLOCATABLE :: CPDISTREF(:)
622:       DOUBLE PRECISION, ALLOCATABLE :: MIN1REDO(:), MIN2REDO(:)622:       DOUBLE PRECISION, ALLOCATABLE :: MIN1REDO(:), MIN2REDO(:)
623:       DOUBLE PRECISION, ALLOCATABLE :: CONDISTREF(:), CONDISTREFLOCAL(:), CONDISTREFON(:), CONDISTREFLOCALON(:)623:       DOUBLE PRECISION, ALLOCATABLE :: CONDISTREF(:), CONDISTREFLOCAL(:), CONDISTREFON(:), CONDISTREFLOCALON(:)
624:       DOUBLE PRECISION, ALLOCATABLE :: CONCUT(:), CONCUTLOCAL(:)624:       DOUBLE PRECISION, ALLOCATABLE :: CONCUT(:), CONCUTLOCAL(:)
625:       DOUBLE PRECISION, ALLOCATABLE :: CONDISTREFFIX(:), CONCUTFIX(:)625:       DOUBLE PRECISION, ALLOCATABLE :: CONDISTREFFIX(:), CONCUTFIX(:)
626:       DOUBLE PRECISION, ALLOCATABLE :: CPCONDISTREF(:)626:       DOUBLE PRECISION, ALLOCATABLE :: CPCONDISTREF(:)
627:       DOUBLE PRECISION, ALLOCATABLE :: CONGEOM(:,:)627:       DOUBLE PRECISION, ALLOCATABLE :: CONGEOM(:,:)
628:       DOUBLE PRECISION, ALLOCATABLE :: LJADDEPS(:,:)628:       DOUBLE PRECISION, ALLOCATABLE :: LJADDEPS(:,:)
629:       DOUBLE PRECISION, ALLOCATABLE :: PYADDEPS(:,:) 
630:       LOGICAL, ALLOCATABLE :: CONACTIVE(:)629:       LOGICAL, ALLOCATABLE :: CONACTIVE(:)
631:       LOGICAL, ALLOCATABLE :: ATOMACTIVE(:)630:       LOGICAL, ALLOCATABLE :: ATOMACTIVE(:)
632:       INTEGER, ALLOCATABLE :: ORDERI(:), ORDERJ(:), REPPOW(:)631:       INTEGER, ALLOCATABLE :: ORDERI(:), ORDERJ(:), REPPOW(:)
633:       INTEGER, ALLOCATABLE :: CONI(:), CONJ(:), CONION(:), CONJON(:)632:       INTEGER, ALLOCATABLE :: CONI(:), CONJ(:), CONION(:), CONJON(:)
634:       INTEGER, ALLOCATABLE :: CONIFIX(:), CONJFIX(:), REPIFIX(:), REPJFIX(:)633:       INTEGER, ALLOCATABLE :: CONIFIX(:), CONJFIX(:), REPIFIX(:), REPJFIX(:)
635:       INTEGER, ALLOCATABLE :: REPI(:), REPJ(:)634:       INTEGER, ALLOCATABLE :: REPI(:), REPJ(:)
636:       INTEGER, ALLOCATABLE :: CPCONI(:), CPCONJ(:)635:       INTEGER, ALLOCATABLE :: CPCONI(:), CPCONJ(:)
637:       INTEGER, ALLOCATABLE :: CPREPI(:), CPREPJ(:)636:       INTEGER, ALLOCATABLE :: CPREPI(:), CPREPJ(:)
638:       DOUBLE PRECISION, ALLOCATABLE :: REPCUT(:), NREPCUT(:), CPREPCUT(:), REPCUTFIX(:)637:       DOUBLE PRECISION, ALLOCATABLE :: REPCUT(:), NREPCUT(:), CPREPCUT(:), REPCUTFIX(:)
639:       INTEGER, ALLOCATABLE :: NREPI(:), NREPJ(:)638:       INTEGER, ALLOCATABLE :: NREPI(:), NREPJ(:)


r31876/gay-berne.f90 2017-02-03 06:30:31.659237676 +0000 r31875/gay-berne.f90 2017-02-03 06:30:32.911254522 +0000
769:                 J2=3*J1769:                 J2=3*J1
770:                 IF(FROZEN(J1)) THEN770:                 IF(FROZEN(J1)) THEN
771:                         G(J2-2)=0.0D0771:                         G(J2-2)=0.0D0
772:                         G(J2-1)=0.0D0772:                         G(J2-1)=0.0D0
773:                         G(J2)=0.0D0773:                         G(J2)=0.0D0
774:                 END IF774:                 END IF
775:         END DO775:         END DO
776:         RETURN776:         RETURN
777:       END SUBROUTINE PYGPERIODIC 777:       END SUBROUTINE PYGPERIODIC 
778: 778: 
779:  
780:  
781: ! PY potential, dc430's implementation 
782: ! with PBC and continuous cutoff added 
783: ! plus extra LJ site 
784: ! 
785: ! This version is addressable with interactions scaled by the input epsilon matrix PYADDEPS 
786: ! 
787:       SUBROUTINE PYGPERIODICADD (X, G, ENERGY, GTEST) 
788:  
789:        use commons, only: BOXLX,BOXLY,BOXLZ,PARAMONOVPBCX,PARAMONOVPBCY,PARAMONOVPBCZ,PCUTOFF,PARAMONOVCUTOFF,& 
790:                 &         natoms,pysignot,pyepsnot,radift,LJSITE,&  
791:                 &         MAXINTERACTIONS,VT,& 
792:                 &         FROZEN, PYADDEPS, NADDTARGET 
793: !                &         PSIGMAATTR, PEPSILONATTR, MYUNIT, PYBINARYTYPE1, PYBINARYT,& 
794: !                &         PSCALEFAC2, PSCALEFAC1, PEPSILON1, BLJSITE, PYA2BIN, PYA1BIN 
795:  
796:        use pymodule 
797:       IMPLICIT NONE 
798:  
799:       DOUBLE PRECISION :: X(3*NATOMS), G(3*NATOMS) 
800:       DOUBLE PRECISION :: RI(3), RJ(3), RIJ(3), NR(3), RIJSQ, ABSRIJ, P(3) 
801:       DOUBLE PRECISION :: AE1(3,3), BE1(3,3), AE2(3,3), BE2(3,3), APB(3,3), APBINV(3,3) 
802:       DOUBLE PRECISION :: FCNT1, FCNT2, SRTFI1, SRTFI2, SRTFI(2), FMIN, LAMDAC1, LAMDAC2, ENERGY 
803:       DOUBLE PRECISION :: RHO1, RHO1SQ, RHO16, RHO112, RHO2, RHO2SQ, RHO26 
804:       DOUBLE PRECISION :: FCTR1, FCTR2, DVDF1, DVDF2  
805:       DOUBLE PRECISION :: DF1PI1, DF1PI2, DF1PI3, DF2PI1, DF2PI2, DF2PI3 
806:       DOUBLE PRECISION :: DF1PJ1, DF1PJ2, DF1PJ3, DF2PJ1, DF2PJ2, DF2PJ3  
807:       DOUBLE PRECISION :: RMI(3,3), RMJ(3,3) 
808:       DOUBLE PRECISION :: DPI1RM(3,3), DPI2RM(3,3), DPI3RM(3,3), DPJ1RM(3,3), DPJ2RM(3,3), DPJ3RM(3,3) 
809:       DOUBLE PRECISION :: DF1DR(3), DF2DR(3), DG1DR(3), DG2DR(3) 
810:       DOUBLE PRECISION :: ARIBRJ(3), XC(3), XCMRI(3), XCMRJ(3), FIJ(3), TIJ(3), TJI(3) 
811:       double precision :: D1ABEZ(3,3), D2ABEZ(3,3), D3ABEZ(3,3), D1ABE(3,3), D2ABE(3,3), D3ABE(3,3)  
812:       LOGICAL          :: GTEST 
813: !     DOUBLE PRECISION :: E(3,3), ESQ(3,3), DE1(3,3), DE2(3,3), DE3(3,3), THETA, THETA2, CT, ST, RM(3,3) 
814:  
815: ! sf344 additions 
816:  
817:       DOUBLE PRECISION :: CLJ1(2),CLJ2(2),CLJ3(2),CLJ4(2),CLJ5(2),CLJ6(2),CLJ7(2),CLJ8(2),CLJ11(2),CLJ12(2), & 
818:                           & CLJ13(2),CLJ14(2),DFDR(2,3),DFP(2,6),dr(6),dCLJ1(2,12),dVDUMMY(12),LJ1(2),DUMMY,DUMMY1,DUMMY2,& 
819:                           & dLJ1(2,12),VDUMMY 
820:       DOUBLE PRECISION :: term2(maxinteractions),term3(maxinteractions), & 
821:                           & xlj(maxinteractions,2,3),rljvec(maxinteractions,3),rljunitvec(maxinteractions,3),& 
822:                           & drlj(maxinteractions,12),rlj(maxinteractions),rlj2(maxinteractions) 
823:       DOUBLE PRECISION :: LLJ(12,maxinteractions), dLLJ1(maxinteractions,12) 
824:       INTEGER          :: k, MJ1, MJ2 
825:  
826:      VT(1:NATOMS/2)=0.0D0 
827:   
828:      term2(:)=1.0D0 
829:      term3(:)=0.0D0 
830:        
831:          ENERGY = 0.D0 
832:          G(:)   = 0.D0 
833:          
834:         DO J1=1, REALNATOMS 
835:             J3      = 3*J1 
836:             J5      = OFFSET + J3 
837:             RI      = X(J3-2:J3) 
838:             P       = X(J5-2:J5) 
839:             angle=sqrt(dot_product(P,P)) 
840:             if(angle>twopi) then 
841: ! normalise angle-axis coordinates 
842:                 X(J5-2:J5)=X(J5-2:J5)/angle 
843:                 do 
844:                   angle=angle-twopi 
845:                   if(angle<2*pi) exit 
846:                 end do 
847: ! multiply with new angle 
848:                 X(J5-2:J5)=X(J5-2:J5)*angle 
849:             end if 
850:  
851:             CALL RMDRVT(P, RMIvec(J1,:,:), DPI1RMvec(J1,:,:), DPI2RMvec(J1,:,:), DPI3RMvec(J1,:,:), GTEST) 
852:  
853:         END DO 
854:  
855:          DO J1 = 1, REALNATOMS 
856:            
857:             MJ1=MOD(J1-1,NADDTARGET)+1 
858:             J3      = 3*J1 
859:             J5      = OFFSET + J3 
860:             RI      = X(J3-2:J3) 
861:             P       = X(J5-2:J5) 
862: !     ROTATION MATRIX 
863:  
864: !            CALL RMDRVT(P, RMI, DPI1RM, DPI2RM, DPI3RM, GTEST) 
865:             RMI(:,:)=RMIvec(J1,:,:) 
866:             DPI1RM(:,:)=DPI1RMvec(J1,:,:) 
867:             DPI2RM(:,:)=DPI2RMvec(J1,:,:) 
868:             DPI3RM(:,:)=DPI3RMvec(J1,:,:) 
869:  
870:             AE1 = MATMUL(RMI,(MATMUL(AEZR1(J1,:,:),(TRANSPOSE(RMI))))) 
871:  
872:             IF (RADIFT) THEN 
873:  
874:                AE2 = MATMUL(RMI,(MATMUL(AEZR2(J1,:,:),(TRANSPOSE(RMI)))))          
875:  
876:             ENDIF 
877:  
878: !     BEGIN INNER LOOP OVER PARTICLES 
879:  
880:             DO J2 = J1 + 1, REALNATOMS 
881:  
882:                MJ2=MOD(J2-1,NADDTARGET)+1 
883:  
884:                J4     = 3*J2 
885:                J6     = OFFSET + J4 
886:                RJ     = X(J4-2:J4)  
887:                P      = X(J6-2:J6) 
888:  
889: !     ROTATION MATRIX 
890:  
891: !               CALL RMDRVT(P, RMJ, DPJ1RM, DPJ2RM, DPJ3RM, GTEST)                
892:                RMJ(:,:)=RMIvec(J2,:,:) 
893:                DPJ1RM(:,:)=DPI1RMvec(J2,:,:) 
894:                DPJ2RM(:,:)=DPI2RMvec(J2,:,:) 
895:                DPJ3RM(:,:)=DPI3RMvec(J2,:,:) 
896:       
897:                BE1 = MATMUL(RMJ,(MATMUL(AEZR1(J2,:,:),(TRANSPOSE(RMJ))))) 
898:  
899:                IF (RADIFT) THEN 
900:     
901:                   BE2 = MATMUL(RMJ,(MATMUL(AEZR2(J2,:,:),(TRANSPOSE(RMJ))))) 
902:  
903:                ENDIF 
904:  
905: !     CALCULATE SEPARATION 
906:  
907:                RIJ    = RI - RJ 
908:                RIJSQ  = DOT_PRODUCT(RIJ,RIJ) 
909:                ABSRIJ = DSQRT(RIJSQ) 
910:                NR     = RIJ / ABSRIJ 
911:  
912:                 IF(PARAMONOVCUTOFF.AND.RIJSQ>cut) GOTO 124 
913:  
914: !     CALCULATE ECF 
915:  
916:                CALL BRENTMIN (0.D0, 0.51D0, 1.D0, AE1, BE1, RIJ, LAMDAC1, FMIN) 
917:  
918:                FCNT1   = - FMIN 
919:                SRTFI1  = 1.D0 / DSQRT(FCNT1) 
920:                APB     = LAMDAC1 * AE1 + (1.D0 - LAMDAC1) * BE1 
921:                 
922:                CALL MTRXIN (APB, APBINV) 
923:  
924:                ARIBRJ =  LAMDAC1 * MATMUL(AE1,RI) + (1.D0 - LAMDAC1) * MATMUL(BE1,RJ) 
925:                XC     =  MATMUL(APBINV, ARIBRJ) 
926:                XCMRI  = XC - RI 
927:                XCMRJ  = XC - RJ 
928:                DF1DR  = - 2.D0 * LAMDAC1 * MATMUL(AE1,XCMRI) 
929:  
930:                D1ABEZ = MATMUL(DPI1RM,AEZR1(J1,:,:)) 
931:                D1ABE  = MATMUL(D1ABEZ,(TRANSPOSE(RMI))) + MATMUL(RMI,(TRANSPOSE(D1ABEZ))) 
932:  
933:                D2ABEZ = MATMUL(DPI2RM,AEZR1(J1,:,:)) 
934:                D2ABE  = MATMUL(D2ABEZ,(TRANSPOSE(RMI))) + MATMUL(RMI,(TRANSPOSE(D2ABEZ))) 
935:  
936:                D3ABEZ = MATMUL(DPI3RM,AEZR1(J1,:,:)) 
937:                D3ABE  = MATMUL(D3ABEZ,(TRANSPOSE(RMI))) + MATMUL(RMI,(TRANSPOSE(D3ABEZ))) 
938:  
939:                DF1PI1 = LAMDAC1*DOT_PRODUCT(XCMRI,MATMUL(D1ABE,XCMRI)) 
940:                DF1PI2 = LAMDAC1*DOT_PRODUCT(XCMRI,MATMUL(D2ABE,XCMRI)) 
941:                DF1PI3 = LAMDAC1*DOT_PRODUCT(XCMRI,MATMUL(D3ABE,XCMRI)) 
942:  
943:                D1ABEZ = MATMUL(DPJ1RM,AEZR1(J2,:,:)) 
944:                D1ABE  = MATMUL(D1ABEZ,(TRANSPOSE(RMJ))) + MATMUL(RMJ,(TRANSPOSE(D1ABEZ))) 
945:  
946:                D2ABEZ = MATMUL(DPJ2RM,AEZR1(J2,:,:)) 
947:                D2ABE  = MATMUL(D2ABEZ,(TRANSPOSE(RMJ))) + MATMUL(RMJ,(TRANSPOSE(D2ABEZ))) 
948:  
949:                D3ABEZ = MATMUL(DPJ3RM,AEZR1(J2,:,:)) 
950:                D3ABE  = MATMUL(D3ABEZ,(TRANSPOSE(RMJ))) + MATMUL(RMJ,(TRANSPOSE(D3ABEZ)))  
951:                 
952:                DF1PJ1 = (1.D0-LAMDAC1)*DOT_PRODUCT(XCMRJ,MATMUL(D1ABE,XCMRJ)) 
953:                DF1PJ2 = (1.D0-LAMDAC1)*DOT_PRODUCT(XCMRJ,MATMUL(D2ABE,XCMRJ)) 
954:                DF1PJ3 = (1.D0-LAMDAC1)*DOT_PRODUCT(XCMRJ,MATMUL(D3ABE,XCMRJ)) 
955:  
956:                RHO1   = PYSIGNOT / (ABSRIJ - ABSRIJ*SRTFI1 + PYSIGNOT) 
957:                RHO1SQ = RHO1*RHO1 
958:                RHO16  = RHO1SQ*RHO1SQ*RHO1SQ 
959:                RHO112 = RHO16 * RHO16 
960:  
961:                FCTR1  = 0.5D0*ABSRIJ*SRTFI1/(FCNT1*PYSIGNOT) 
962:                DG1DR  = (1.D0-SRTFI1)*NR/PYSIGNOT + FCTR1*DF1DR 
963:                DVDF1  = -2.D0*RHO112*RHO1*FCTR1 
964:  
965:                IF (RADIFT) THEN 
966:  
967:                   CALL BRENTMIN (0.D0, 0.51D0, 1.D0, AE2, BE2, RIJ, LAMDAC2, FMIN) 
968:  
969:                   FCNT2   = - FMIN 
970:                   SRTFI2  = 1.D0 / DSQRT(FCNT2) 
971:                   APB     = LAMDAC2 * AE2 + (1.D0 - LAMDAC2) * BE2 
972:  
973:                   CALL MTRXIN (APB, APBINV) 
974:  
975:                   ARIBRJ =  LAMDAC2 * MATMUL(AE2,RI) + (1.D0 - LAMDAC2) * MATMUL(BE2,RJ) 
976:                   XC     =  MATMUL(APBINV, ARIBRJ) 
977:                   XCMRI  = XC - RI 
978:                   XCMRJ  = XC - RJ 
979:                   DF2DR  = - 2.D0 * LAMDAC2 * MATMUL(AE2,XCMRI) 
980:  
981:                   RHO2   = PYSIGNOT / (ABSRIJ - ABSRIJ*SRTFI2 + PYSIGNOT) 
982:                   RHO2SQ = RHO2*RHO2 
983:                   RHO26  = RHO2SQ*RHO2SQ*RHO2SQ 
984:                 
985:                   FCTR2  = 0.5D0*ABSRIJ*SRTFI2/(FCNT2*PYSIGNOT) 
986:                   DG2DR  = (1.D0-SRTFI2)*NR/PYSIGNOT+FCTR2*DF2DR 
987:                   DVDF2  = RHO26*RHO2*FCTR2 
988:  
989:                   D1ABEZ = MATMUL(DPI1RM,AEZR2(J1,:,:)) 
990:                   D1ABE  = MATMUL(D1ABEZ,(TRANSPOSE(RMI))) + MATMUL(RMI,(TRANSPOSE(D1ABEZ))) 
991:  
992:                   D2ABEZ = MATMUL(DPI2RM,AEZR2(J1,:,:)) 
993:                   D2ABE  = MATMUL(D2ABEZ,(TRANSPOSE(RMI))) + MATMUL(RMI,(TRANSPOSE(D2ABEZ))) 
994:  
995:                   D3ABEZ = MATMUL(DPI3RM,AEZR2(J1,:,:)) 
996:                   D3ABE  = MATMUL(D3ABEZ,(TRANSPOSE(RMI))) + MATMUL(RMI,(TRANSPOSE(D3ABEZ))) 
997:  
998:                   DF2PI1 = LAMDAC2*DOT_PRODUCT(MATMUL(XCMRI,D1ABE),XCMRI) 
999:                   DF2PI2 = LAMDAC2*DOT_PRODUCT(MATMUL(XCMRI,D2ABE),XCMRI) 
1000:                   DF2PI3 = LAMDAC2*DOT_PRODUCT(MATMUL(XCMRI,D3ABE),XCMRI) 
1001:  
1002:                   D1ABEZ = MATMUL(DPJ1RM,AEZR2(J2,:,:)) 
1003:                   D1ABE  = MATMUL(D1ABEZ,(TRANSPOSE(RMJ))) + MATMUL(RMJ,(TRANSPOSE(D1ABEZ))) 
1004:  
1005:                   D2ABEZ = MATMUL(DPJ2RM,AEZR2(J2,:,:)) 
1006:                   D2ABE  = MATMUL(D2ABEZ,(TRANSPOSE(RMJ))) + MATMUL(RMJ,(TRANSPOSE(D2ABEZ))) 
1007:  
1008:                   D3ABEZ = MATMUL(DPJ3RM,AEZR2(J2,:,:)) 
1009:                   D3ABE  = MATMUL(D3ABEZ,(TRANSPOSE(RMJ))) + MATMUL(RMJ,(TRANSPOSE(D3ABEZ))) 
1010:  
1011:                   DF2PJ1 = (1.D0-LAMDAC2)*DOT_PRODUCT(MATMUL(XCMRJ,D1ABE),XCMRJ) 
1012:                   DF2PJ2 = (1.D0-LAMDAC2)*DOT_PRODUCT(MATMUL(XCMRJ,D2ABE),XCMRJ) 
1013:                   DF2PJ3 = (1.D0-LAMDAC2)*DOT_PRODUCT(MATMUL(XCMRJ,D3ABE),XCMRJ) 
1014:  
1015:                ELSE 
1016:  
1017:                   RHO2   = RHO1 
1018:                   RHO26  = RHO16 
1019:                   DG2DR  = DG1DR 
1020:                   DVDF2  = RHO26*RHO2*FCTR1 
1021:                   DF2PI1 = DF1PI1 
1022:                   DF2PI2 = DF1PI2 
1023:                   DF2PI3 = DF1PI3 
1024:                   DF2PJ1 = DF1PJ1 
1025:                   DF2PJ2 = DF1PJ2 
1026:                   DF2PJ3 = DF1PJ3 
1027:  
1028:                   SRTFI2 = SRTFI1 
1029:                   DF2DR  = DF1DR 
1030:                   RHO2SQ = RHO1SQ 
1031:  
1032:                ENDIF 
1033:  
1034:  
1035: ! Correction terms to the potential if we require a cutoff at rc.  
1036:                  SRTFI(1)=SRTFI1 
1037:                  SRTFI(2)=SRTFI2 
1038:                  DFP(1,1)=DF1PI1 
1039:                  DFP(1,2)=DF1PI2 
1040:                  DFP(1,3)=DF1PI3 
1041:                  DFP(1,4)=DF1PJ1 
1042:                  DFP(1,5)=DF1PJ2 
1043:                  DFP(1,6)=DF1PJ3 
1044:                  DFP(2,1)=DF2PI1 
1045:                  DFP(2,2)=DF2PI2 
1046:                  DFP(2,3)=DF2PI3 
1047:                  DFP(2,4)=DF2PJ1 
1048:                  DFP(2,5)=DF2PJ2 
1049:                  DFP(2,6)=DF2PJ3 
1050:                  DFDR(1,:)=DF1DR(:) 
1051:                  DFDR(2,:)=DF2DR(:) 
1052:                  LJ1(1)=RHO1 
1053:                  LJ1(2)=RHO2 
1054:  !     CALCULATE PY POTENTIAL ENERGY 
1055: 123 CONTINUE  
1056:         VDUMMY = 0.0D0 
1057:                   
1058:                 VDUMMY=0.0D0 
1059:                IF(LJSITE) THEN   
1060:                 do k=1,maxinteractions ! k=1 -- interaction between repulsive primary 'apex' sites 
1061:                          ! k=2 and k=3 -- interaction between secondary and primary 'apex' sites 
1062:                          ! k=4 -- interaction between secondary 'apex' sites (normal LJ interaction)  
1063:                         ! trying to modify code to allow for binary systems.  
1064:                         ! apex site heights will be defined in absolute units, 
1065:                         ! hence PYA1bin(J1,1) etc. will be removed from below 
1066:  
1067:                    IF(k==1) THEN 
1068:                         DUMMY1=PSCALEFAC1vec(J1)!*PYA1bin(J1,1) 
1069:                         DUMMY2=PSCALEFAC1vec(J2)!*PYA1bin(J2,1) 
1070:                    ELSE IF(k==2) THEN 
1071:                         DUMMY1=PSCALEFAC1vec(J1)!*PYA1bin(J1,1) 
1072:                         DUMMY2=-PSCALEFAC2vec(J2)!*PYA1bin(J2,1) 
1073:                    ELSE IF(k==3) THEN 
1074:                         DUMMY1=-PSCALEFAC2vec(J1)!*PYA1bin(J1,1) 
1075:                         DUMMY2=PSCALEFAC1vec(J2)!*PYA1bin(J2,1) 
1076:                    ELSE 
1077:                         DUMMY1=-PSCALEFAC2vec(J1)!*PYA1bin(J1,1) 
1078:                         DUMMY2=-PSCALEFAC2vec(J2)!*PYA1bin(J2,1) 
1079:                    END IF 
1080:                         ! first particle 
1081:                         xlj(k,1,:)=RI+DUMMY1*MATMUL(RMI,vecsbf)    ! vecsbf: (1,0,0) in the body frame of ellipsoid 
1082:  
1083:                         ! second particle 
1084:                         xlj(k,2,:)=RJ+DUMMY2*MATMUL(RMJ,vecsbf) 
1085:  
1086:                         ! separation between the LJ sites 
1087:                         rlj2(k)=(xlj(k,2,1)-xlj(k,1,1))**2+(xlj(k,2,2)-xlj(k,1,2))**2+(xlj(k,2,3)-xlj(k,1,3))**2 
1088:                         rlj(k)=sqrt(rlj2(k)) 
1089:                         rljvec(k,1)=xlj(k,2,1)-xlj(k,1,1) 
1090:                         rljvec(k,2)=xlj(k,2,2)-xlj(k,1,2) 
1091:                         rljvec(k,3)=xlj(k,2,3)-xlj(k,1,3) 
1092:  
1093:                         DUMMY=1.0D0/rlj(k) 
1094:                         rljunitvec(k,:)=rljvec(k,:)*DUMMY !/rlj(k) 
1095:  
1096:                         drlj(k,1)=DUMMY*(xlj(k,2,1)-xlj(k,1,1))         !drlj/dx1 
1097:                         drlj(k,2)=DUMMY*(xlj(k,2,2)-xlj(k,1,2))         !drlj/dy1 
1098:                         drlj(k,3)=DUMMY*(xlj(k,2,3)-xlj(k,1,3))         !drlj/dz1 
1099:                         drlj(k,4)=-drlj(k,1)                               !drlj/dx2 
1100:                         drlj(k,5)=-drlj(k,2)                               !drlj/dy2 
1101:                         drlj(k,6)=-drlj(k,3)                               !drlj/dz2 
1102:                         drlj(k,7) =-DUMMY*DUMMY1*DOT_PRODUCT(rljvec(k,:),MATMUL(DPI1RM,vecsbf)) !drlj/dpx1 
1103:                         drlj(k,8) =-DUMMY*DUMMY1*DOT_PRODUCT(rljvec(k,:),MATMUL(DPI2RM,vecsbf)) !drlj/dpy1 
1104:                         drlj(k,9) =-DUMMY*DUMMY1*DOT_PRODUCT(rljvec(k,:),MATMUL(DPI3RM,vecsbf)) !drlj/dpz1 
1105:                         drlj(k,10) =  DUMMY*DUMMY2*DOT_PRODUCT(rljvec(k,:),MATMUL(DPJ1RM,vecsbf)) !drlj/dpx2 
1106:                         drlj(k,11) =  DUMMY*DUMMY2*DOT_PRODUCT(rljvec(k,:),MATMUL(DPJ2RM,vecsbf)) !drlj/dpy2 
1107:                         drlj(k,12) =  DUMMY*DUMMY2*DOT_PRODUCT(rljvec(k,:),MATMUL(DPJ3RM,vecsbf)) !drlj/dpz2 
1108:  
1109:               ! interaction between the extra LJ sites: 
1110:                         LLJ(1,k)=sigma1(k)*DUMMY !/rlj(k) 
1111:                         LLJ(2,k)=LLJ(1,k)**2 
1112:                         LLJ(3,k)=LLJ(2,k)*LLJ(1,k) 
1113:                         LLJ(4,k)=LLJ(2,k)**2 
1114:                         LLJ(5,k)=LLJ(4,k)*LLJ(1,k) 
1115:                         LLJ(6,k)=LLJ(4,k)*LLJ(2,k) 
1116:                         LLJ(7,k)=LLJ(6,k)*LLJ(1,k) 
1117:                         LLJ(11,k)=LLJ(5,k)*LLJ(6,k) 
1118:                         LLJ(12,k)=LLJ(6,k)*LLJ(6,k) 
1119:  
1120:                             DO j=1,12 
1121:                                 dLLJ1(k,j) =-sigma1(k)*DUMMY*DUMMY*drlj(k,j) 
1122:                             END DO 
1123:  
1124:                  VDUMMY=VDUMMY+4.0D0*epsilon1(k,J1,J2)*term2(k)*(LLJ(12,k) - attr(k)*LLJ(6,k)) ! extra LJ sites (12-6) 
1125:               end do ! k=1,4 
1126:              END IF ! IF(LJSITE) 
1127:  
1128:                 VDUMMY = VDUMMY + 4.0D0 * PYEPSNOT * (RHO112 - 1.0D0 * RHO26)*PYADDEPS(MJ2,MJ1) 
1129:  
1130:                ENERGY = ENERGY + VDUMMY 
1131:                VT(J1) = VT(J1) + VDUMMY 
1132:                VT(J2) = VT(J2) + VDUMMY        ! pair potentials 
1133:  
1134:  
1135:             dVDUMMY(:) = 0.0D0 
1136:  
1137: !     CALCULATE GRADIENTPYADDEPS(MJ2,MJ1) 
1138:  
1139:              FIJ = 0.0D0 
1140:              TIJ = 0.0D0 
1141:              TJI = 0.0D0 
1142:  
1143:              FIJ = 0.0D0 
1144:              TIJ = 0.0D0 
1145:              TJI = 0.0D0 
1146:              dVDUMMY(:)=0.0D0 
1147:  
1148:            IF(LJSITE) THEN 
1149:             do k=1,maxinteractions 
1150:              do j=1,3 
1151:                dVDUMMY(j) = dVDUMMY(j) + 4.0D0*epsilon1(k,J1,J2)*(12.0D0*LLJ(11,k)*dLLJ1(k,j))*term2(k) + & 
1152:                           & 4.0D0*epsilon1(k,J1,J2)*LLJ(12,k)*term3(k)*drlj(k,j)! extra LJ site derivatives, currently 
1153:                                                                                 ! only repulsive (LJ12) 
1154:                dVDUMMY(j) = dVDUMMY(j) - attr(k)*(4.0D0*epsilon1(k,J1,J2)*(6.0D0*LLJ(5,k)*dLLJ1(k,j))*term2(k) + & 
1155:                           & 4.0D0*epsilon1(k,J1,J2)*LLJ(6,k)*term3(k)*drlj(k,j)) ! attractive secondary apex site 
1156:                FIJ(j) = dVDUMMY(j) 
1157:              end do 
1158:              do j=7,12 
1159:                dVDUMMY(j) = dVDUMMY(j) + 4.0D0*epsilon1(k,J1,J2)*(12.0D0*LLJ(11,k)*dLLJ1(k,j))*term2(k) + & 
1160:                           & 4.0D0*epsilon1(k,J1,J2)*LLJ(12,k)*term3(k)*drlj(k,j)! extra LJ site derivatives, currently 
1161:                                                                                 ! only repulsive 
1162:                dVDUMMY(j) = dVDUMMY(j) - attr(k)*(4.0D0*epsilon1(k,J1,J2)*(6.0D0*LLJ(5,k)*dLLJ1(k,j))*term2(k) + & 
1163:                           & 4.0D0*epsilon1(k,J1,J2)*LLJ(6,k)*term3(k)*drlj(k,j)) ! attractive secondary apex site 
1164:              end do 
1165:             end do 
1166:            END IF !if LJSITE  
1167:                FIJ    = FIJ + 24.0D0*PYEPSNOT*(2.D0*RHO112*RHO1*DG1DR - 1.0D0 * RHO26*RHO2*DG2DR)*PYADDEPS(MJ2,MJ1) 
1168:                TIJ(1) = DVDF1*DF1PI1 + 1.0D0* DVDF2*DF2PI1 
1169:                TIJ(2) = DVDF1*DF1PI2 + 1.0D0*  DVDF2*DF2PI2 
1170:                TIJ(3) = DVDF1*DF1PI3 + 1.0D0*  DVDF2*DF2PI3 
1171:                TJI(1) = DVDF1*DF1PJ1 + 1.0D0*  DVDF2*DF2PJ1 
1172:                TJI(2) = DVDF1*DF1PJ2 + 1.0D0*  DVDF2*DF2PJ2 
1173:                TJI(3) = DVDF1*DF1PJ3 + 1.0D0*  DVDF2*DF2PJ3 
1174:                TIJ = dVDUMMY(7:9) + 24.0D0 * PYEPSNOT * TIJ*PYADDEPS(MJ2,MJ1) 
1175:                TJI = dVDUMMY(10:12) + 24.0D0 * PYEPSNOT * TJI*PYADDEPS(MJ2,MJ1) 
1176:  
1177:                G(J3-2:J3) = G(J3-2:J3) - FIJ 
1178:                G(J4-2:J4) = G(J4-2:J4) + FIJ 
1179:                G(J5-2:J5) = G(J5-2:J5) + TIJ 
1180:                G(J6-2:J6) = G(J6-2:J6) + TJI 
1181:  
1182: !     END INNER LOOP OVER PARTICLES 
1183: 124 CONTINUE 
1184:             ENDDO 
1185:  
1186: !     END OUTER LOOP OVER PARTICLES 
1187:  
1188:          ENDDO 
1189:         DO J1=1,NATOMS 
1190:                 J2=3*J1 
1191:                 IF(FROZEN(J1)) THEN 
1192:                         G(J2-2)=0.0D0 
1193:                         G(J2-1)=0.0D0 
1194:                         G(J2)=0.0D0 
1195:                 END IF 
1196:         END DO 
1197:         RETURN 
1198:       END SUBROUTINE PYGPERIODICADD 
1199:  
1200:  
1201: SUBROUTINE ECF(OVERLAPTEST,ECFvalue,x1,x2,y1,y2,z1,z2,px11,px12,py11,py12,pz11,pz12,a,b,c)779: SUBROUTINE ECF(OVERLAPTEST,ECFvalue,x1,x2,y1,y2,z1,z2,px11,px12,py11,py12,pz11,pz12,a,b,c)
1202: 780: 
1203: 781: 
1204: !Calculates the value of the ellipsoid contact function782: !Calculates the value of the ellipsoid contact function
1205: 783: 
1206: !use commons784: !use commons
1207: IMPLICIT NONE785: IMPLICIT NONE
1208: 786: 
1209: DOUBLE PRECISION, intent(in)       :: x1,x2,y1,y2,z1,z2,px11,px12,py11,py12,pz11,pz12,a,b,c787: DOUBLE PRECISION, intent(in)       :: x1,x2,y1,y2,z1,z2,px11,px12,py11,py12,pz11,pz12,a,b,c
1210: DOUBLE PRECISION, intent(out)      :: ECFvalue788: DOUBLE PRECISION, intent(out)      :: ECFvalue


r31876/io1.f 2017-02-03 06:30:31.899240905 +0000 r31875/io1.f 2017-02-03 06:30:33.231258829 +0000
629:       ELSE IF (WATERKZT) THEN629:       ELSE IF (WATERKZT) THEN
630: 630: 
631:          WRITE(MYUNIT,'(I4,A)') NATOMS/2,' polarisableKZ water molecules '631:          WRITE(MYUNIT,'(I4,A)') NATOMS/2,' polarisableKZ water molecules '
632: 632: 
633:       ELSE IF (GAYBERNET) THEN633:       ELSE IF (GAYBERNET) THEN
634: 634: 
635:          WRITE(MYUNIT,'(I4,A)') NATOMS/2,' Gay-Berne atoms '635:          WRITE(MYUNIT,'(I4,A)') NATOMS/2,' Gay-Berne atoms '
636: 636: 
637:       ELSE IF (PYGPERIODICT) THEN637:       ELSE IF (PYGPERIODICT) THEN
638: 638: 
639: 639:          WRITE(MYUNIT,'(I4,A)') NATOMS/2,' ellipsoids '
640:          IF  (PYADDT) THEN 
641:             IF (SORTT) THEN 
642:                WRITE(MYUNIT,'(A)') 'Turning off SORT option for PYADD' 
643:                SORTT=.FALSE. 
644:             ENDIF 
645:             WRITE(MYUNIT,'(I4,A)') NATOMS/2,' addressable ellipsoids' 
646:          ELSE 
647:             WRITE(MYUNIT,'(I4,A)') NATOMS/2,' ellipsoids ' 
648:          ENDIF 
649: 640: 
650:       ELSE IF (PYT) THEN641:       ELSE IF (PYT) THEN
651: 642: 
652:          WRITE(MYUNIT,'(I4,A)') NATOMS/2,' PY building blocks'643:          WRITE(MYUNIT,'(I4,A)') NATOMS/2,' PY building blocks'
653: 644: 
654:       ELSE IF (GAYBERNEDCT) THEN645:       ELSE IF (GAYBERNEDCT) THEN
655: 646: 
656:          WRITE(MYUNIT,'(I4,A)') NATOMS/2,' Gay-Berne atoms '647:          WRITE(MYUNIT,'(I4,A)') NATOMS/2,' Gay-Berne atoms '
657:          DO J2=1,NPAR648:          DO J2=1,NPAR
658:             WRITE(MYUNIT,*) COORDS(:,NPAR)649:             WRITE(MYUNIT,*) COORDS(:,NPAR)


r31876/keywords.f 2017-02-03 06:30:32.135244081 +0000 r31875/keywords.f 2017-02-03 06:30:33.455261843 +0000
1161: !1161: !
1162:       MLQT=.FALSE.1162:       MLQT=.FALSE.
1163:       MLQPROB=.FALSE.1163:       MLQPROB=.FALSE.
1164:       MLQDONE=.FALSE.1164:       MLQDONE=.FALSE.
1165:       MLQNORM=.FALSE.1165:       MLQNORM=.FALSE.
1166:       MLQLAMBDA=0.0D01166:       MLQLAMBDA=0.0D0
1167:       MLQSTART=11167:       MLQSTART=1
1168: 1168: 
1169:       LJADDT=.FALSE.1169:       LJADDT=.FALSE.
1170:       LJADD2T=.FALSE.1170:       LJADD2T=.FALSE.
1171:       PYADDT=.FALSE. 
1172: 1171: 
1173:       DUMPMQT=.FALSE.1172:       DUMPMQT=.FALSE.
1174:       1173:       
1175:       CALL FILE_OPEN('data', DATA_UNIT, .FALSE.)1174:       CALL FILE_OPEN('data', DATA_UNIT, .FALSE.)
1176:       1175:       
1177: !      OPEN (5,FILE='data',STATUS='OLD')1176: !      OPEN (5,FILE='data',STATUS='OLD')
1178: 1177: 
1179: !190   CALL INPUT(END,5)1178: !190   CALL INPUT(END,5)
1180: 190   CALL INPUT(END, DATA_UNIT)1179: 190   CALL INPUT(END, DATA_UNIT)
1181:       IF (.NOT. END) THEN1180:       IF (.NOT. END) THEN
4434:          JUMPMOVE(IX)=.TRUE.4433:          JUMPMOVE(IX)=.TRUE.
4435:          CALL READI(JUMPTO(IX))4434:          CALL READI(JUMPTO(IX))
4436:          JDUMP(JUMPTO(IX))=.TRUE.4435:          JDUMP(JUMPTO(IX))=.TRUE.
4437:          IF (NITEMS.GT.3) CALL READI(JUMPINT(IX))4436:          IF (NITEMS.GT.3) CALL READI(JUMPINT(IX))
4438: 4437: 
4439:       ELSE IF (WORD.EQ.'LB2') THEN4438:       ELSE IF (WORD.EQ.'LB2') THEN
4440:          LB2T=.TRUE.4439:          LB2T=.TRUE.
4441: !4440: !
4442: !  Addressable LJ4441: !  Addressable LJ
4443: !4442: !
4444:       ELSE IF (WORD.EQ.'LJADD') THEN4443:          ELSE IF (WORD.EQ.'LJADD') THEN
4445:          LJADDT=.TRUE.4444:             LJADDT=.TRUE.
4446:          LUNIT=GETUNIT()4445:             LUNIT=GETUNIT()
4447:          OPEN(LUNIT,FILE='epsilon',STATUS='OLD')4446:             OPEN(LUNIT,FILE='epsilon',STATUS='OLD')
4448:          IF (.NOT.ALLOCATED(LJADDEPS)) ALLOCATE(LJADDEPS(NATOMS,NATOMS))4447:             IF (.NOT.ALLOCATED(LJADDEPS)) ALLOCATE(LJADDEPS(NATOMS,NATOMS))
4449:          DO J1=1,NATOMS4448:             DO J1=1,NATOMS
4450:             DO J2=1,NATOMS4449:                DO J2=1,NATOMS
4451:                READ(LUNIT,*) LJADDEPS(J2,J1)4450:                   READ(LUNIT,*) LJADDEPS(J2,J1)
4452:                WRITE(MYUNIT,'(2I6,G20.10)') J1,J2,LJADDEPS(J2,J1)4451:                   WRITE(MYUNIT,'(2I6,G20.10)') J1,J2,LJADDEPS(J2,J1)
 4452:                ENDDO
4453:             ENDDO4453:             ENDDO
4454:          ENDDO4454:             CLOSE(LUNIT)
4455:          CLOSE(LUNIT)4455:          ELSE IF (WORD.EQ.'LJADD2') THEN
4456:       ELSE IF (WORD.EQ.'LJADD2') THEN4456:             LJADDT=.TRUE.
4457:          LJADDT=.TRUE.4457:             LJADD2T=.TRUE.
4458:          LJADD2T=.TRUE.4458:             CALL READI(NADDTARGET)
4459:          CALL READI(NADDTARGET)4459:             WRITE(MYUNIT,'(A,I6)') 'keyword> Target cluster size is ',NADDTARGET
4460:          WRITE(MYUNIT,'(A,I6)') 'keyword> Target cluster size is ',NADDTARGET4460:             IF (MOD(NATOMS,NADDTARGET).NE.0) THEN
4461:          IF (MOD(NATOMS,NADDTARGET).NE.0) THEN4461:                WRITE(MYUNIT,'(A,I6)') 'keyword> ERROR, target cluster size is not a factor of the nummber of the atoms ',NATOMS
4462:             WRITE(MYUNIT,'(A,I6)') 'keyword> ERROR, target cluster size is not a factor of the number of the atoms ',NATOMS4462:                STOP
4463:             STOP4463:             ENDIF
4464:          ENDIF4464:             LUNIT=GETUNIT()
4465:          LUNIT=GETUNIT()4465:             OPEN(LUNIT,FILE='epsilon',STATUS='OLD')
4466:          OPEN(LUNIT,FILE='epsilon',STATUS='OLD')4466:             IF (.NOT.ALLOCATED(LJADDEPS)) ALLOCATE(LJADDEPS(NADDTARGET,NADDTARGET))
4467:          IF (.NOT.ALLOCATED(LJADDEPS)) ALLOCATE(LJADDEPS(NADDTARGET,NADDTARGET))4467:             DO J1=1,NADDTARGET
4468:          DO J1=1,NADDTARGET4468:                DO J2=1,NADDTARGET
4469:             DO J2=1,NADDTARGET4469:                   READ(LUNIT,*) LJADDEPS(J2,J1)
4470:                READ(LUNIT,*) LJADDEPS(J2,J1)4470:                   WRITE(MYUNIT,'(2I6,G20.10)') J1,J2,LJADDEPS(J2,J1)
4471:                WRITE(MYUNIT,'(2I6,G20.10)') J1,J2,LJADDEPS(J2,J1)4471:                ENDDO
4472:             ENDDO4472:             ENDDO
4473:          ENDDO4473:             CLOSE(LUNIT)
4474:          CLOSE(LUNIT) 
4475: !4474: !
4476: ! LJAT4475: ! LJAT
4477: !4476: !
4478:       ELSE IF (WORD.EQ.'LJAT') THEN4477:       ELSE IF (WORD.EQ.'LJAT') THEN
4479:          LJATT=.TRUE.4478:          LJATT=.TRUE.
4480:          CALL READF(ZSTAR)4479:          CALL READF(ZSTAR)
4481:          LJATTOC=2.423D04480:          LJATTOC=2.423D0
4482:          IF (NITEMS.EQ.3) CALL READF(LJATTOC)4481:          IF (NITEMS.EQ.3) CALL READF(LJATTOC)
4483:       ELSE IF (WORD.EQ.'LOCALSAMPLE') THEN4482:       ELSE IF (WORD.EQ.'LOCALSAMPLE') THEN
4484:          LOCALSAMPLET=.TRUE.4483:          LOCALSAMPLET=.TRUE.
5816:       5815:       
5817:        ELSE IF (WORD.EQ.'ROTATEHINGE') THEN5816:        ELSE IF (WORD.EQ.'ROTATEHINGE') THEN
5818:           ROTATEHINGET=.TRUE.5817:           ROTATEHINGET=.TRUE.
5819: !         Read ROTATEHINGEFREQ 5818: !         Read ROTATEHINGEFREQ 
5820:           IF (NITEMS.GT.1) CALL READI(ROTATEHINGEFREQ)5819:           IF (NITEMS.GT.1) CALL READI(ROTATEHINGEFREQ)
5821: !         Read in ROTATEFACTOR5820: !         Read in ROTATEFACTOR
5822:           IF (NITEMS.GT.2) CALL READF(HINGEROTATEFACTOR)5821:           IF (NITEMS.GT.2) CALL READF(HINGEROTATEFACTOR)
5823:           WRITE(MYUNIT,'(A)') ' keywords> Enabled hinge rotation moves for the plate potential'5822:           WRITE(MYUNIT,'(A)') ' keywords> Enabled hinge rotation moves for the plate potential'
5824:           WRITE(MYUNIT,'(A,I2,A)') ' keywords> Hinges will be rotated every ',ROTATEHINGEFREQ,' steps'5823:           WRITE(MYUNIT,'(A,I2,A)') ' keywords> Hinges will be rotated every ',ROTATEHINGEFREQ,' steps'
5825:           WRITE(MYUNIT,'(A,F20.10)') ' keywords> Hinge angle rotatefactor =',HINGEROTATEFACTOR5824:           WRITE(MYUNIT,'(A,F20.10)') ' keywords> Hinge angle rotatefactor =',HINGEROTATEFACTOR
5826:           CALL HINGE_INITIALISE5825:           CALL HINGE_INITIALISE
5827:          5826:          
5828: !5827: !
5829: ! csw34> Rigid body rotation moves. Each rigid body is randomly rotated about its COM every ROTATERIGIDFREQ steps.5828: ! csw34> Rigid body rotation moves. Each rigid body is randomly rotated about its COM every ROTATERIGIDFREQ steps.
5830: !        ROTATEFACTOR scales the maximum rotation with 1.0 being complete freedom to rotate.5829: !        ROTATEFACTOR scales the maximum rotation with 1.0 being complete freedom to rotate.
5831: !5830: !
5832:       ELSE IF (WORD.EQ.'ROTATERIGID') THEN5831:       ELSE IF (WORD.EQ.'ROTATERIGID') THEN
5833:          ROTATERIGIDT=.TRUE.5832:          ROTATERIGIDT=.TRUE.
5834: ! Read ROTATERIGIDFREQ 5833: ! Read ROTATERIGIDFREQ 
5835:          IF (NITEMS.GT.1) CALL READI(ROTATERIGIDFREQ)5834:          IF (NITEMS.GT.1) CALL READI(ROTATERIGIDFREQ)
5836: ! Read in ROTATEFACTOR5835: ! Read in ROTATEFACTOR
6708:          ELLIPSOIDT=.TRUE.6707:          ELLIPSOIDT=.TRUE.
6709:          RIGID=.TRUE.6708:          RIGID=.TRUE.
6710:          NRBSITES=16709:          NRBSITES=1
6711:          CALL READF(GBANISOTROPYR)6710:          CALL READF(GBANISOTROPYR)
6712:          CALL READF(GBWELLDEPTHR)6711:          CALL READF(GBWELLDEPTHR)
6713:          CALL READF(PSIGMA0(1))6712:          CALL READF(PSIGMA0(1))
6714:          CALL READF(PEPSILON0)6713:          CALL READF(PEPSILON0)
6715:          CALL READF(GBMU)6714:          CALL READF(GBMU)
6716:          CALL READF(GBNU)6715:          CALL READF(GBNU)
6717:          ALLOCATE(SITE(NRBSITES,3))6716:          ALLOCATE(SITE(NRBSITES,3))
6718: ! 
6719: ! Addressable PY 
6720: ! 
6721:       ELSE IF (WORD.EQ.'PYADD') THEN 
6722:          PYADDT=.TRUE. 
6723:          CALL READI(NADDTARGET) 
6724:          WRITE(MYUNIT,'(A,I6)') 'keyword> Target cluster size is ',NADDTARGET 
6725:          IF (MOD(NATOMSALLOC/2,NADDTARGET).NE.0) THEN 
6726:             WRITE(MYUNIT,'(A,I6)') 'keyword> ERROR, target cluster size is not a factor of the number of PY particles ',  
6727:      &          NATOMSALLOC/2   
6728:             STOP 
6729:          ENDIF 
6730:          LUNIT=GETUNIT() 
6731:          OPEN(LUNIT,FILE='epsilon',STATUS='OLD') 
6732:          IF (.NOT.ALLOCATED(PYADDEPS)) ALLOCATE(PYADDEPS(NADDTARGET,NADDTARGET)) 
6733:          DO J1=1,NADDTARGET 
6734:             DO J2=1,NADDTARGET 
6735:                READ(LUNIT,*) PYADDEPS(J2,J1) 
6736:                WRITE(MYUNIT,'(2I6,G20.10)') J1,J2,PYADDEPS(J2,J1) 
6737:             ENDDO 
6738:          ENDDO 
6739:          CLOSE(LUNIT) 
6740: 6717: 
6741:       ELSE IF (WORD .EQ.'PYGPERIODIC') THEN6718:       ELSE IF (WORD .EQ.'PYGPERIODIC') THEN
6742: 6719: 
6743:          PYGPERIODICT = .TRUE.6720:          PYGPERIODICT = .TRUE.
6744:          ELLIPSOIDT = .TRUE.6721:          ELLIPSOIDT = .TRUE.
6745:          RIGID = .TRUE.6722:          RIGID = .TRUE.
6746:          CALL READF(PYA1(1))6723:          CALL READF(PYA1(1))
6747:          CALL READF(PYA1(2))6724:          CALL READF(PYA1(2))
6748:          CALL READF(PYA1(3))6725:          CALL READF(PYA1(3))
6749:          CALL READF(PYA2(1))6726:          CALL READF(PYA2(1))


r31876/potential.f90 2017-02-03 06:30:32.427248010 +0000 r31875/potential.f90 2017-02-03 06:30:33.687264964 +0000
1124:    ELSE IF (WATERDCT) THEN1124:    ELSE IF (WATERDCT) THEN
1125:       CALL WATERPDC (X, GRAD, EREAL, GRADT)1125:       CALL WATERPDC (X, GRAD, EREAL, GRADT)
1126: 1126: 
1127:    ELSE IF (WATERKZT) THEN1127:    ELSE IF (WATERKZT) THEN
1128:       CALL WATERPKZ (X, GRAD, EREAL, GRADT)1128:       CALL WATERPKZ (X, GRAD, EREAL, GRADT)
1129: 1129: 
1130:    ELSE IF (GAYBERNET) THEN1130:    ELSE IF (GAYBERNET) THEN
1131:       CALL GAYBERNE(X, GRAD, EREAL, GRADT,.FALSE.)1131:       CALL GAYBERNE(X, GRAD, EREAL, GRADT,.FALSE.)
1132: 1132: 
1133:    ELSE IF (PYGPERIODICT .OR. PYBINARYT) THEN1133:    ELSE IF (PYGPERIODICT .OR. PYBINARYT) THEN
1134:       IF (PYADDT) THEN1134:       CALL PYGPERIODIC(X, GRAD, EREAL, GRADT)
1135:          CALL PYGPERIODICADD(X, GRAD, EREAL, GRADT) 
1136:       ELSE 
1137:          CALL PYGPERIODIC(X, GRAD, EREAL, GRADT) 
1138:       ENDIF 
1139: 1135: 
1140:    ELSE IF (LJCAPSIDT) THEN1136:    ELSE IF (LJCAPSIDT) THEN
1141:       CALL LJCAPSIDMODEL(X, GRAD, EREAL, GRADT)1137:       CALL LJCAPSIDMODEL(X, GRAD, EREAL, GRADT)
1142: 1138: 
1143:    ELSE IF (STICKYT) THEN1139:    ELSE IF (STICKYT) THEN
1144:       CALL RADR(X, GRAD, EREAL, GRADT)1140:       CALL RADR(X, GRAD, EREAL, GRADT)
1145:       CALL STICKY(X, GRAD, EREAL, GRADT,.FALSE.)1141:       CALL STICKY(X, GRAD, EREAL, GRADT,.FALSE.)
1146: !      DIFF=1.0D-31142: !      DIFF=1.0D-3
1147: !      WRITE(*, *) 'analytic and numerical gradients:'1143: !      WRITE(*, *) 'analytic and numerical gradients:'
1148: !      DO J1=1, 3*NATOMS1144: !      DO J1=1, 3*NATOMS


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0