hdiff output

r30777/SBM.f 2016-07-15 07:30:21.207295771 +0100 r30776/SBM.f 2016-07-15 07:30:21.455299109 +0100
  8:   8: 
  9:       LOGICAL :: CALLED=.FALSE.  9:       LOGICAL :: CALLED=.FALSE.
 10:       LOGICAL GTEST, STEST 10:       LOGICAL GTEST, STEST
 11:  11: 
 12:         if(.NOT.CALLED)then 12:         if(.NOT.CALLED)then
 13:  13: 
 14:         write(*,*) 14:         write(*,*)
 15:         write(*,*)  'Calculations will use a Structure-based SMOG model:' 15:         write(*,*)  'Calculations will use a Structure-based SMOG model:'
 16:         write(*,*)  '   For more information on SMOG models, see:' 16:         write(*,*)  '   For more information on SMOG models, see:'
 17:         write(*,*)  '   Software: Noel, et al. PLoS Comput Biol 12, e1004794, 2016.' 17:         write(*,*)  '   Software: Noel, et al. PLoS Comput Biol 12, e1004794, 2016.'
 18:         write(*,*)  '   Model: see reference list at smog-server.org/refs.html' 18:         write(*,*)  '   Model: Whitford, et al. Prot Struct Func Bioinfo 75, 430-441, 2009.'
 19:         write(*,*) 19:         write(*,*)
 20:  20: 
 21:         call SBMinit(NATOMS) 21:         call SBMinit(NATOMS)
 22:         CALLED=.TRUE. 22:         CALLED=.TRUE.
 23:         SBMHESSATOM=-1 
 24:         endIF 23:         endIF
 25:  24: 
 26:  25: 
 27: ! call the energy routine 26: ! call the energy routine
 28:       call calc_energy_SBM(qo,natoms,GRAD,energy) 27:       call calc_energy_SBM(qo,natoms,GRAD,energy)
 29:  28: 
 30:       IF (STEST) THEN 29:       IF (STEST) THEN
 31:          PRINT '(A)','ERROR - second derivatives not available' 30:          PRINT '(A)','ERROR - second derivatives not available'
 32:          STOP 31:          STOP
 33:       ENDIF 32:       ENDIF
 34:       return 33:       return
 35:       end 34:       end
 36:  35: 
 37: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 36: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 38: !* SBMinit 37: !* SBMinit() reads the atom positions from file.  If 1 is selected for *
  38: !* startt then the velocities are assigned, otherwise, they are read   *
  39: !* by selecting 2, or generated by selecting 3                         *
 39: !*********************************************************************** 40: !***********************************************************************
 40:  41: 
 41:       subroutine SBMINIT(NATOMS) 42:       subroutine SBMINIT(NATOMS)
 42:       USE KEY 43:       USE KEY
 43:       USE SBMDATA 44:       USE SBMDATA
 44:       USE COMMONS, only: ATMASS 45:       USE COMMONS, only: ATMASS
 45:       USE GENRIGID, only: RIGIDINIT 46:       USE GENRIGID, only: RIGIDINIT
 46:       implicit NONE 47:       implicit NONE
 47:  48: 
 48:         integer i,j,K,MaxCon,NATOMS,storage, dummy,ANr, Ib22,Ib21,IT1,JT1, 49:         integer i,j,K,MaxCon,NATOMS,storage, dummy,ANr, Ib22,Ib21,IT1,JT1,
 49:      Q KT1,IT2,JT2,KT2,IP1,JP1,KP1,LP1,IP2,JP2,KP2,LP2,nBA1,nTA1,nPA1,  50:      Q KT1,IT2,JT2,KT2,IP1,JP1,KP1,LP1,IP2,JP2,KP2,LP2,nBA1,nTA1,nPA1, 
 50:      Q nBA2,nTA2,nPA2,ind1,ind2,ANt,MDT1, MDT2, cl1,cl2,tempi, 51:      Q nBA2,nTA2,nPA2,ind1,ind2,ANt,MDT1, MDT2, cl1,cl2,tempi,
 51:      Q ATYPE,POS 52:      Q ATYPE,POS
 52:       DOUBLE PRECISION :: Sigma, EpsC 53:       DOUBLE PRECISION :: Sigma, EpsC
 53:       DOUBLE PRECISION,ALLOCATABLE,DIMENSION(:) :: TMPSTT 54:       DOUBLE PRECISION,ALLOCATABLE,DIMENSION(:) :: TMPSTT
 54:       DOUBLE PRECISION RDIFF,ALPHA 55:       DOUBLE PRECISION RDIFF,ALPHA
 55:       INTEGER BONDTYPE,MAXBONDSPERATOM 56:       INTEGER BONDTYPE
 56:       INTEGER M 57:       INTEGER M
 57:       integer AA,BB,ANTEMP 58:       integer AA,BB,ANTEMP
 58:       DOUBLE PRECISION dx,dy,dz 59:       DOUBLE PRECISION dx,dy,dz
 59:       double precision PI,ARG1,ARG2,ARG3,ARG4 60:       double precision PI,ARG1,ARG2,ARG3,ARG4
 60:       DOUBLE PRECISION RSig, Reps,DC 61:       DOUBLE PRECISION RSig, Reps,DC
 61:       integer TMPARG 62:       integer TMPARG
 62:       character TMPCHAR 63:       character TMPCHAR
 63:       integer TMPINT,nexc,I1,I2 64:       integer TMPINT,nexc,I1,I2
 64:       INTEGER TT1,TT2 
 65:       double precision TMPREAL,concentration 65:       double precision TMPREAL,concentration
 66:       LOGICAL TARR(NATOMS) 66:       LOGICAL TARR(NATOMS)
 67:       INTEGER EXCLUSIONS,NUMOFEXCLUSIONS 67:       INTEGER EXCLUSIONS,NUMOFEXCLUSIONS
 68:       INTEGER MAXEXCLUSIONS,MAXEXCLUSIONSELEC 68:       INTEGER MAXEXCLUSIONS,MAXEXCLUSIONSELEC
 69:       PARAMETER(MAXEXCLUSIONS=20) 69:       PARAMETER(MAXEXCLUSIONS=20)
 70:       PARAMETER(MAXEXCLUSIONSELEC=5) 70:       PARAMETER(MAXEXCLUSIONSELEC=5)
 71:       DIMENSION EXCLUSIONS(NATOMS*MAXEXCLUSIONS),NUMOFEXCLUSIONS(NATOMS) 71:       DIMENSION EXCLUSIONS(NATOMS*MAXEXCLUSIONS),NUMOFEXCLUSIONS(NATOMS)
 72:       INTEGER IOSTATUS,J1,J2,NRIGIDBODY,ATOMID 72:       INTEGER IOSTATUS,J1,J2,NRIGIDBODY,ATOMID
 73:       CHARACTER(LEN=10) CHECK1 73:       CHARACTER(LEN=10) CHECK1
 74: !$    INTEGER NTHREADS,OMP_GET_NUM_THREADS,OMP_GET_THREAD_NUM,TID  74: !$    INTEGER NTHREADS,OMP_GET_NUM_THREADS,OMP_GET_THREAD_NUM,TID 
229: 229: 
230:           read(30,*)230:           read(30,*)
231:           read(30,*) NBA231:           read(30,*) NBA
232:           write(*,*) NBA, 'bonds'232:           write(*,*) NBA, 'bonds'
233:           ALLOCATE(IB1(NBA))233:           ALLOCATE(IB1(NBA))
234:           ALLOCATE(IB2(NBA))234:           ALLOCATE(IB2(NBA))
235:           ALLOCATE(RB(NBA))235:           ALLOCATE(RB(NBA))
236:           ALLOCATE(BK(NBA))236:           ALLOCATE(BK(NBA))
237:         do i=1, nBA237:         do i=1, nBA
238:           read(30,*) Ib1(i), Ib2(i),bondtype,Rb(i), bK(i)238:           read(30,*) Ib1(i), Ib2(i),bondtype,Rb(i), bK(i)
239:  
240:           !!! ADD BONDS OF TYPE 6239:           !!! ADD BONDS OF TYPE 6
241:           IF(max(IB1(I),IB2(I))-min(IB1(I),IB2(I)) .gt. MAXSEP .AND. bondtype .EQ. 1)THEN240:           IF(max(IB1(I),IB2(I))-min(IB1(I),IB2(I)) .gt. MAXSEP .AND. bondtype .EQ. 1)THEN
242:             WRITE(*,*) 'WARNING: DIFFERENCE IN BONDED ATOM INDICES IS GREATER THAN MAXSEP'241:             WRITE(*,*) 'WARNING: DIFFERENCE IN BONDED ATOM INDICES IS GREATER THAN MAXSEP'
243:             WRITE(*,*) '         THIS CAN LEAD TO POOR PERFORMANCE.'242:             WRITE(*,*) '         THIS CAN LEAD TO POOR PERFORMANCE.'
244:             WRITE(*,*) '         CONSIDER RECOMPILING WITH A LARGER VALUE OF MAXSEP.'243:             WRITE(*,*) '         CONSIDER RECOMPILING WITH A LARGER VALUE OF MAXSEP.'
245:           ENDIF244:           ENDIF
246:           IF(bondtype .eq. 1 .or. bondtype .eq. 6)THEN245:           IF(bondtype .eq. 1 .or. bondtype .eq. 6)THEN
247:             CALL INCLUDEEXCLUSIONS(NATOMS,IB1(I),IB2(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,.TRUE.)          246:             CALL INCLUDEEXCLUSIONS(NATOMS,IB1(I),IB2(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,.TRUE.)          
248:           ELSEIF(bondtype .eq. 6)THEN247:           ELSEIF(bondtype .eq. 6)THEN
249:             CALL INCLUDEEXCLUSIONS(NATOMS,IB1(I),IB2(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,.FALSE.)          248:             CALL INCLUDEEXCLUSIONS(NATOMS,IB1(I),IB2(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,.FALSE.)          
250:           ENDIF249:           ENDIF
251:         end do250:         end do
252: 251: 
253: ! organizing for rapid hessian evaluations 
254:         ALLOCATE(SBMNB(NATOMS)) 
255:         do i=1,NATOMS  
256:           SBMNB(I)=0 
257:         enddo 
258:  
259:         do i=1, nBA 
260:           SBMNB(IB1(I))=SBMNB(IB1(I))+1 
261:           SBMNB(IB2(I))=SBMNB(IB2(I))+1 
262:         enddo 
263:         MAXBONDSPERATOM=0 
264:         do i=1,NATOMS 
265:           IF(SBMNB(I) .GT. MAXBONDSPERATOM)THEN 
266:                 MAXBONDSPERATOM=SBMNB(I) 
267:           ENDIF 
268:         ENDDO 
269:           ALLOCATE(SBMBONDLIST(NATOMS,MAXBONDSPERATOM)) 
270:  
271:         do i=1,NATOMS  
272:           SBMNB(I)=0 
273:         enddo 
274:  
275:         do i=1, nBA 
276:           TT1=IB1(I) 
277:           TT2=IB2(I) 
278:           SBMNB(TT1)=SBMNB(TT1)+1 
279:           SBMNB(TT2)=SBMNB(TT2)+1 
280:           SBMBONDLIST(TT1,SBMNB(TT1))=I 
281:           SBMBONDLIST(TT2,SBMNB(TT2))=I 
282:         enddo 
283:  
284:           read(30,*)252:           read(30,*)
285:           read(30,*) NTA253:           read(30,*) NTA
286:           write(*,*) NTA, 'angles'254:           write(*,*) NTA, 'angles'
287:         ALLOCATE(IT(NTA))255:         ALLOCATE(IT(NTA))
288:         ALLOCATE(JT(NTA))256:         ALLOCATE(JT(NTA))
289:         ALLOCATE(KT(NTA))257:         ALLOCATE(KT(NTA))
290:         ALLOCATE(ANTC(NTA))258:         ALLOCATE(ANTC(NTA))
291:         ALLOCATE(TK(NTA))259:         ALLOCATE(TK(NTA))
292:         do i=1, nTA260:         do i=1, nTA
293:           read(30,*) IT(i), JT(i), KT(i), ANTC(i), TK(i)261:           read(30,*) IT(i), JT(i), KT(i), ANTC(i), TK(i)
529: 497: 
530:       subroutine SBMBonds(grad,energy, natoms)498:       subroutine SBMBonds(grad,energy, natoms)
531:       USE KEY499:       USE KEY
532:       USE SBMDATA500:       USE SBMDATA
533:       implicit NONE501:       implicit NONE
534:       integer I2,J2,I23,J23,I,N,J,NATOMS502:       integer I2,J2,I23,J23,I,N,J,NATOMS
535:       DOUBLE PRECISION grad(3*NATOMS),energy503:       DOUBLE PRECISION grad(3*NATOMS),energy
536:       DOUBLE PRECISION r2, f, r1504:       DOUBLE PRECISION r2, f, r1
537:       DOUBLE PRECISION dx,dy,dz505:       DOUBLE PRECISION dx,dy,dz
538: 506: 
539:       IF(SBMHESSATOM .EQ. -1)THEN  
540: !$OMP PARALLEL PRIVATE(I,I2,J2,I23,J23,DX,DY,DZ,R1,R2,F)REDUCTION(+:ENERGY,grad)507: !$OMP PARALLEL PRIVATE(I,I2,J2,I23,J23,DX,DY,DZ,R1,R2,F)REDUCTION(+:ENERGY,grad)
541: !$OMP DO508: !$OMP DO
542:         DO I=1, nBA509:         DO I=1, nBA
543:            I2 = Ib1(I)510:            I2 = Ib1(I)
544:            J2 = Ib2(I)511:            J2 = Ib2(I)
545:              I23=I2*3512:            I23=I2*3
546:              J23=J2*3513:            J23=J2*3
547:              dx = XSBM(I2) - XSBM(J2)514:            dx = XSBM(I2) - XSBM(J2)
548:              dy = YSBM(I2) - YSBM(J2)515:            dy = YSBM(I2) - YSBM(J2)
549:              dz = ZSBM(I2) - ZSBM(J2)516:            dz = ZSBM(I2) - ZSBM(J2)
550:              r2 = dx**2 + dy**2 + dz**2517:            r2 = dx**2 + dy**2 + dz**2
551:              r1 = sqrt(r2)518:            r1 = sqrt(r2)
552:              ENERGY = ENERGY + bk(I)*(r1-Rb(I))**2/2.0519:            ENERGY = ENERGY + bk(I)*(r1-Rb(I))**2/2.0
553:              f = -bk(I)*(r1-Rb(I))/r1520:            f = -bk(I)*(r1-Rb(I))/r1
554:              grad(I23-2) = grad(I23-2) - f * dx521:            grad(I23-2) = grad(I23-2) - f * dx
555:              grad(I23-1) = grad(I23-1) - f * dy522:            grad(I23-1) = grad(I23-1) - f * dy
556:              grad(I23)   = grad(I23)   - f * dz523:            grad(I23)   = grad(I23)   - f * dz
557:              grad(J23-2) = grad(J23-2) + f * dx524:            grad(J23-2) = grad(J23-2) + f * dx
558:              grad(J23-1) = grad(J23-1) + f * dy525:            grad(J23-1) = grad(J23-1) + f * dy
559:              grad(J23)   = grad(J23)   + f * dz526:            grad(J23)   = grad(J23)   + f * dz
560:        ENDDO 527:        ENDDO 
561: !$OMP END DO528: !$OMP END DO
562: !$OMP END PARALLEL529: !$OMP END PARALLEL
563: 530: 
564:       ELSE 
565:         DO I=1, SBMNB(SBMHESSATOM) 
566:            I2 = IB1(SBMBONDLIST(SBMHESSATOM,I)) 
567:            J2 = IB2(SBMBONDLIST(SBMHESSATOM,I)) 
568:              I23=I2*3 
569:              J23=J2*3 
570:              dx = XSBM(I2) - XSBM(J2) 
571:              dy = YSBM(I2) - YSBM(J2) 
572:              dz = ZSBM(I2) - ZSBM(J2) 
573:              r2 = dx**2 + dy**2 + dz**2 
574:              r1 = sqrt(r2) 
575:              ENERGY = ENERGY + bk(I)*(r1-Rb(I))**2/2.0 
576:              f = -bk(I)*(r1-Rb(I))/r1 
577:              grad(I23-2) = grad(I23-2) - f * dx 
578:              grad(I23-1) = grad(I23-1) - f * dy 
579:              grad(I23)   = grad(I23)   - f * dz 
580:              grad(J23-2) = grad(J23-2) + f * dx 
581:              grad(J23-1) = grad(J23-1) + f * dy 
582:              grad(J23)   = grad(J23)   + f * dz 
583:        ENDDO  
584:  
585:       ENDIF 
586:  
587:       END531:       END
588: 532: 
589: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^END OF SBMBONDS^^^^^^^^^^^^^^^^^^^^^^^^^^^^^533: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^END OF SBMBONDS^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
590: 534: 
591: 535: 
592: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<536: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
593: !* SBMANGL  computes the Force due to the bond angles                   *537: !* SBMANGL  computes the Force due to the bond angles                   *
594: !* This code is modeled after how AMBER performs angle forces           *538: !* This code is modeled after how AMBER performs angle forces           *
595: !***********************************************************************539: !***********************************************************************
596: 540: 
613:         ONE=1.0557:         ONE=1.0
614:         NEGONE=-1.0558:         NEGONE=-1.0
615: !$OMP PARALLEL PRIVATE(I3,J3,K3,I33,J33,K33,RIJ,RKJ,RIK,ANT,XIJ,YIJ,ZIJ,XKJ,YKJ,559: !$OMP PARALLEL PRIVATE(I3,J3,K3,I33,J33,K33,RIJ,RKJ,RIK,ANT,XIJ,YIJ,ZIJ,XKJ,YKJ,
616: !$OMP&  ZKJ,DF,CT0,CT1,CT2,DA,ST,CIK,CII,CKK,DT1,DT2,DT3,DT4,560: !$OMP&  ZKJ,DF,CT0,CT1,CT2,DA,ST,CIK,CII,CKK,DT1,DT2,DT3,DT4,
617: !$OMP&  DT5,DT6,DT7,DT8,DT9,STH) REDUCTION(+:ENERGY,grad)561: !$OMP&  DT5,DT6,DT7,DT8,DT9,STH) REDUCTION(+:ENERGY,grad)
618: !$OMP DO562: !$OMP DO
619:           DO II = 1, nTA563:           DO II = 1, nTA
620:             I3 = IT(II)564:             I3 = IT(II)
621:             J3 = JT(II)565:             J3 = JT(II)
622:             K3 = KT(II)566:             K3 = KT(II)
623:             IF(SBMHESSATOM .EQ. -1 .OR. SBMHESSATOM .EQ. I3 .OR.567:             I33=I3*3
624:      Q        SBMHESSATOM .EQ. J3 .OR. SBMHESSATOM .EQ. K3)THEN568:             J33=J3*3
625:               I33=I3*3569:             K33=K3*3
626:               J33=J3*3570:             XIJ = XSBM(I3)-XSBM(J3)
627:               K33=K3*3571:             YIJ = YSBM(I3)-YSBM(J3)
628:               XIJ = XSBM(I3)-XSBM(J3)572:             ZIJ = ZSBM(I3)-ZSBM(J3)
629:               YIJ = YSBM(I3)-YSBM(J3)573:             XKJ = XSBM(K3)-XSBM(J3)
630:               ZIJ = ZSBM(I3)-ZSBM(J3)574:             YKJ = YSBM(K3)-YSBM(J3)
631:               XKJ = XSBM(K3)-XSBM(J3)575:             ZKJ = ZSBM(K3)-ZSBM(J3)
632:               YKJ = YSBM(K3)-YSBM(J3)576:             RIJ = XIJ*XIJ+YIJ*YIJ+ZIJ*ZIJ
633:               ZKJ = ZSBM(K3)-ZSBM(J3)577:             RKJ = XKJ*XKJ+YKJ*YKJ+ZKJ*ZKJ
634:               RIJ = XIJ*XIJ+YIJ*YIJ+ZIJ*ZIJ578:             RIK = SQRT(RIJ*RKJ)
635:               RKJ = XKJ*XKJ+YKJ*YKJ+ZKJ*ZKJ579:             CT0 = (XIJ*XKJ+YIJ*YKJ+ZIJ*ZKJ)/RIK
636:               RIK = SQRT(RIJ*RKJ)580:             CT1 = MAX(NEGONE,CT0)
637:               CT0 = (XIJ*XKJ+YIJ*YKJ+ZIJ*ZKJ)/RIK581:             CT2 = MIN(ONE,CT1)
638:               CT1 = MAX(NEGONE,CT0)582:             ANT = ACOS(CT2)
639:               CT2 = MIN(ONE,CT1)583:             DA = ANT - ANTC(II)
640:               ANT = ACOS(CT2)584:             DF = TK(II)*DA
641:               DA = ANT - ANTC(II)585:             ST = -(DF)/SIN(ANT)
642:               DF = TK(II)*DA586:             STH = ST*CT2
643:               ST = -(DF)/SIN(ANT)587:             CIK = ST/RIK
644:               STH = ST*CT2588:             CII = STH/RIJ
645:               CIK = ST/RIK589:             CKK = STH/RKJ
646:               CII = STH/RIJ590:             DT1 = CIK*XKJ-CII*XIJ
647:               CKK = STH/RKJ591:             DT2 = CIK*YKJ-CII*YIJ
648:               DT1 = CIK*XKJ-CII*XIJ592:             DT3 = CIK*ZKJ-CII*ZIJ
649:               DT2 = CIK*YKJ-CII*YIJ593:             DT7 = CIK*XIJ-CKK*XKJ
650:               DT3 = CIK*ZKJ-CII*ZIJ594:             DT8 = CIK*YIJ-CKK*YKJ
651:               DT7 = CIK*XIJ-CKK*XKJ595:             DT9 = CIK*ZIJ-CKK*ZKJ
652:               DT8 = CIK*YIJ-CKK*YKJ596:             DT4 = -DT1-DT7
653:               DT9 = CIK*ZIJ-CKK*ZKJ597:             DT5 = -DT2-DT8
654:               DT4 = -DT1-DT7598:             DT6 = -DT3-DT9
655:               DT5 = -DT2-DT8599:             grad(I33-2) = grad(I33-2)+ DT1
656:               DT6 = -DT3-DT9600:             grad(I33-1) = grad(I33-1)+ DT2
657:               grad(I33-2) = grad(I33-2)+ DT1601:             grad(I33)   = grad(I33)  + DT3
658:               grad(I33-1) = grad(I33-1)+ DT2602:             grad(J33-2) = grad(J33-2)+ DT4
659:               grad(I33)   = grad(I33)  + DT3603:             grad(J33-1) = grad(J33-1)+ DT5
660:               grad(J33-2) = grad(J33-2)+ DT4604:             grad(J33)   = grad(J33)  + DT6
661:               grad(J33-1) = grad(J33-1)+ DT5605:             grad(K33-2) = grad(K33-2)+ DT7
662:               grad(J33)   = grad(J33)  + DT6606:             grad(K33-1) = grad(K33-1)+ DT8
663:               grad(K33-2) = grad(K33-2)+ DT7607:             grad(K33)   = grad(K33)  + DT9
664:               grad(K33-1) = grad(K33-1)+ DT8608:             ENERGY = ENERGY + TK(II)*(ANTC(II)- ANT)**2/2.0
665:               grad(K33)   = grad(K33)  + DT9 
666:               ENERGY = ENERGY + TK(II)*(ANTC(II)- ANT)**2/2.0 
667:             ENDIF 
668:           END DO609:           END DO
669: !$OMP END DO610: !$OMP END DO
670: !$OMP END PARALLEL611: !$OMP END PARALLEL
671: 612: 
672:        RETURN613:        RETURN
673:        END614:        END
674: 615: 
675: !^^^^^^^^^^^^^^^^^^^^^^^^End of SBMANGL^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^616: !^^^^^^^^^^^^^^^^^^^^^^^^End of SBMANGL^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
676: 617: 
677: 618: 
720: !$OMP& S,HGoverG,FGoverG,A1,A3,TT1,TT2,TT3,TT4,TT1X,TT1Y,TT1Z,TT2X,661: !$OMP& S,HGoverG,FGoverG,A1,A3,TT1,TT2,TT3,TT4,TT1X,TT1Y,TT1Z,TT2X,
721: !$OMP& TT2Y,TT2Z,TT3X,TT3Y,TT3Z,TT4X,TT4Y,TT4Z) REDUCTION(+:energy,grad)662: !$OMP& TT2Y,TT2Z,TT3X,TT3Y,TT3Z,TT4X,TT4Y,TT4Z) REDUCTION(+:energy,grad)
722: !!663: !!
723: !$OMP DO664: !$OMP DO
724:           DO II = 1,nPA665:           DO II = 1,nPA
725: 666: 
726:             I3 = IP(II)667:             I3 = IP(II)
727:             J3 = JP(II)668:             J3 = JP(II)
728:             K3 = KP(II)669:             K3 = KP(II)
729:             L3 = LP(II)670:             L3 = LP(II)
730:             IF(SBMHESSATOM .EQ. -1 .OR. SBMHESSATOM .EQ. I3 .OR.671:             I33=I3*3
731:      Q        SBMHESSATOM .EQ. J3 .OR. SBMHESSATOM .EQ. K3 .OR. 672:             J33=J3*3
732:      Q        SBMHESSATOM .EQ. L3)THEN 673:             K33=K3*3
733:               I33=I3*3674:             L33=L3*3
734:               J33=J3*3675:             XIJ = XSBM(I3)-XSBM(J3)
735:               K33=K3*3676:             YIJ = YSBM(I3)-YSBM(J3)
736:               L33=L3*3677:             ZIJ = ZSBM(I3)-ZSBM(J3)
737:               XIJ = XSBM(I3)-XSBM(J3)678:             XKJ = XSBM(K3)-XSBM(J3)
738:               YIJ = YSBM(I3)-YSBM(J3)679:             YKJ = YSBM(K3)-YSBM(J3)
739:               ZIJ = ZSBM(I3)-ZSBM(J3)680:             ZKJ = ZSBM(K3)-ZSBM(J3)
740:               XKJ = XSBM(K3)-XSBM(J3)681:             RKJ = sqrt(XKJ**2+YKJ**2+ZKJ**2)
741:               YKJ = YSBM(K3)-YSBM(J3)682:             XKL = XSBM(K3)-XSBM(L3)
742:               ZKJ = ZSBM(K3)-ZSBM(J3)683:             YKL = YSBM(K3)-YSBM(L3)
743:               RKJ = sqrt(XKJ**2+YKJ**2+ZKJ**2)684:             ZKL = ZSBM(K3)-ZSBM(L3)                                  
744:               XKL = XSBM(K3)-XSBM(L3)685:             FGoverG=-(XIJ*XKJ+YIJ*YKJ+ZIJ*ZKJ)/RKJ
745:               YKL = YSBM(K3)-YSBM(L3)686:             HGoverG=(XKL*XKJ+YKL*YKJ+ ZKL*ZKJ)/RKJ
746:               ZKL = ZSBM(K3)-ZSBM(L3)                                   
747:               FGoverG=-(XIJ*XKJ+YIJ*YKJ+ZIJ*ZKJ)/RKJ 
748:               HGoverG=(XKL*XKJ+YKL*YKJ+ ZKL*ZKJ)/RKJ 
749: C DX is the M vector and G is the N vector687: C DX is the M vector and G is the N vector
750:               DX = YIJ*ZKJ-ZIJ*YKJ688:             DX = YIJ*ZKJ-ZIJ*YKJ
751:               DY = ZIJ*XKJ-XIJ*ZKJ689:             DY = ZIJ*XKJ-XIJ*ZKJ
752:               DZ = XIJ*YKJ-YIJ*XKJ690:             DZ = XIJ*YKJ-YIJ*XKJ
753:               GX = ZKJ*YKL-YKJ*ZKL691:             GX = ZKJ*YKL-YKJ*ZKL
754:               GY = XKJ*ZKL-ZKJ*XKL692:             GY = XKJ*ZKL-ZKJ*XKL
755:               GZ = YKJ*XKL-XKJ*YKL693:             GZ = YKJ*XKL-XKJ*YKL
756:               FXI = SQRT(DX*DX+DY*DY+DZ*DZ)694:             FXI = SQRT(DX*DX+DY*DY+DZ*DZ)
757:               FYI = SQRT(GX*GX+GY*GY+GZ*GZ)695:             FYI = SQRT(GX*GX+GY*GY+GZ*GZ)
758:               CT = DX*GX+DY*GY+DZ*GZ696:             CT = DX*GX+DY*GY+DZ*GZ
759:               z10 = 1.0/FXI697:             z10 = 1.0/FXI
760:               z20 = 1.0/FYI698:             z20 = 1.0/FYI
761:               Z12 = Z10*Z20699:             Z12 = Z10*Z20
762:               Z1 = Z10700:             Z1 = Z10
763:               Z2 = Z20701:             Z2 = Z20
764:               CT0 = MIN(ONE,CT*Z12)702:             CT0 = MIN(ONE,CT*Z12)
765:               CT1 = MAX(NEGONE,CT0)703:             CT1 = MAX(NEGONE,CT0)
766:               S = XKJ*(DZ*GY-DY*GZ)+YKJ*(DX*GZ-DZ*GX)+ZKJ*(DY*GX-DX*GY)704:             S = XKJ*(DZ*GY-DY*GZ)+YKJ*(DX*GZ-DZ*GX)+ZKJ*(DY*GX-DX*GY)
767:               AP0 = ACOS(CT1)705:             AP0 = ACOS(CT1)
768:               AP1 = PI-SIGN(AP0,S)706:             AP1 = PI-SIGN(AP0,S)
769:               CT = AP1707:             CT = AP1
770:               CPHI = COS(AP1)708:             CPHI = COS(AP1)
771:               SPHI = SIN(AP1)709:             SPHI = SIN(AP1)
772: ! Here is the energy part710: ! Here is the energy part
773:             A1=CT-PHISBM(II)711:           A1=CT-PHISBM(II)
774:             A3=A1*3712:           A3=A1*3
775:   713: 
776:           if(PHITYPE(II) .eq. 1)then714:         if(PHITYPE(II) .eq. 1)then
777:             energy =  energy + PK(II)*(3.0/2.0-cos(A1)-0.5*cos(A3))715:           energy =  energy + PK(II)*(3.0/2.0-cos(A1)-0.5*cos(A3))
778: ! dE/dPHI 
779:               DF=PK(II)*(sin(A1)+1.5*sin(A3)) 
780:           elseif(PHITYPE(II) .eq. 2)then 
781:             if(A1 .gt. PI)then 
782:                   A1=A1-2*PI 
783:             elseif(A1 .lt. -PI)then 
784:                   A1=A1+2*PI 
785:             endif 
786:    
787:             energy =  energy + 0.5*PK(II)*A1**2 
788: ! dE/dPHI716: ! dE/dPHI
789:               DF=PK(II)*A1717:             DF=PK(II)*(sin(A1)+1.5*sin(A3))
790:   718:         elseif(PHITYPE(II) .eq. 2)then
791:           else719:           if(A1 .gt. PI)then
792:   720:                 A1=A1-2*PI
793:           write(*,*) 'unrecognized dihedral type', PHITYPE(II)721:           elseif(A1 .lt. -PI)then
794:           STOP722:                 A1=A1+2*PI
795:           endif723:           endif
796:   724: 
 725:           energy =  energy + 0.5*PK(II)*A1**2
 726: ! dE/dPHI
 727:             DF=PK(II)*A1
 728: 
 729:         else
 730: 
 731:         write(*,*) 'unrecognized dihedral type', PHITYPE(II)
 732:         STOP
 733:         endif
 734: 
797: ! insert the new 735: ! insert the new 
798:   736: 
799: ! now, do dPhi/dX737: ! now, do dPhi/dX
800:   738: 
801: ! |G|/|A|**2 739: ! |G|/|A|**2 
802:             TT1 = Z1*Z1*RKJ*DF740:           TT1 = Z1*Z1*RKJ*DF
803: ! FG/(A**2*|G|)741: ! FG/(A**2*|G|)
804:             TT2 = FGoverG*Z1*Z1*DF742:           TT2 = FGoverG*Z1*Z1*DF
805: ! HG/(B**2*|G|)743: ! HG/(B**2*|G|)
806:             TT3 = HGoverG*Z2*Z2*DF744:           TT3 = HGoverG*Z2*Z2*DF
807: ! |G|/|B|**2 745: ! |G|/|B|**2 
808:             TT4 = Z2*Z2*RKJ*DF746:           TT4 = Z2*Z2*RKJ*DF
809: ! note: negatives are flipped from paper because A=-DX747: ! note: negatives are flipped from paper because A=-DX
810:             TT1X=TT1*DX748:           TT1X=TT1*DX
811:             TT1Y=TT1*DY749:           TT1Y=TT1*DY
812:             TT1Z=TT1*DZ750:           TT1Z=TT1*DZ
813:             TT2X=TT2*DX751:           TT2X=TT2*DX
814:             TT2Y=TT2*DY752:           TT2Y=TT2*DY
815:             TT2Z=TT2*DZ753:           TT2Z=TT2*DZ
816:             TT3X=TT3*GX754:           TT3X=TT3*GX
817:             TT3Y=TT3*GY755:           TT3Y=TT3*GY
818:             TT3Z=TT3*GZ756:           TT3Z=TT3*GZ
819:             TT4X=TT4*GX757:           TT4X=TT4*GX
820:             TT4Y=TT4*GY758:           TT4Y=TT4*GY
821:             TT4Z=TT4*GZ759:           TT4Z=TT4*GZ
822:             grad(I33-2) =  grad(I33-2)  + TT1X  760:           grad(I33-2) =  grad(I33-2)  + TT1X  
823:             grad(I33-1) =  grad(I33-1)  + TT1Y761:           grad(I33-1) =  grad(I33-1)  + TT1Y
824:             grad(I33)   =  grad(I33)    + TT1Z762:           grad(I33)   =  grad(I33)    + TT1Z
825:             grad(J33-2) =  grad(J33-2)  - TT1X - TT2X - TT3X763:           grad(J33-2) =  grad(J33-2)  - TT1X - TT2X - TT3X
826:             grad(J33-1) =  grad(J33-1)  - TT1Y - TT2Y - TT3Y764:           grad(J33-1) =  grad(J33-1)  - TT1Y - TT2Y - TT3Y
827:             grad(J33)   =  grad(J33)    - TT1Z - TT2Z - TT3Z765:           grad(J33)   =  grad(J33)    - TT1Z - TT2Z - TT3Z
828:             grad(K33-2) =  grad(K33-2)  + TT2X + TT3X - TT4X766:           grad(K33-2) =  grad(K33-2)  + TT2X + TT3X - TT4X
829:             grad(K33-1) =  grad(K33-1)  + TT2Y + TT3Y - TT4Y767:           grad(K33-1) =  grad(K33-1)  + TT2Y + TT3Y - TT4Y
830:             grad(K33)   =  grad(K33)    + TT2Z + TT3Z - TT4Z768:           grad(K33)   =  grad(K33)    + TT2Z + TT3Z - TT4Z
831:             grad(L33-2) =  grad(L33-2)  + TT4X769:           grad(L33-2) =  grad(L33-2)  + TT4X
832:             grad(L33-1) =  grad(L33-1)  + TT4Y770:           grad(L33-1) =  grad(L33-1)  + TT4Y
833:             grad(L33)   =  grad(L33)    + TT4Z771:           grad(L33)   =  grad(L33)    + TT4Z
834:           ENDIF 
835:         END DO772:         END DO
836: !$OMP END DO 773: !$OMP END DO 
837: !$OMP END PARALLEL774: !$OMP END PARALLEL
838: 775: 
839:           END776:           END
840: 777: 
841: 778: 
842: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^END of SBMDihedral^^^^^^^^^^^^^^^^^^^^^^^^^^^779: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^END of SBMDihedral^^^^^^^^^^^^^^^^^^^^^^^^^^^
843: 780: 
844: 781: 
864: ! type 2 is 12-12801: ! type 2 is 12-12
865: ! type 5 as gaussian, no ex-vol802: ! type 5 as gaussian, no ex-vol
866: ! type 6 as gaussian, w ex-vol803: ! type 6 as gaussian, w ex-vol
867: ! type 7 as dual gaussian804: ! type 7 as dual gaussian
868: !$OMP PARALLEL PRIVATE(I,CO1,CO2,CO3,CO4,CO5,CO6,C1,C2,C13,C23,DX,DY,DZ,805: !$OMP PARALLEL PRIVATE(I,CO1,CO2,CO3,CO4,CO5,CO6,C1,C2,C13,C23,DX,DY,DZ,
869: !$OMP& R2,RM2,RM6,RM10,RM12,R1,E1,G1,G2,E1D,G1D,G2D,F_OVER_R,ET) REDUCTION(+:ENERGY,grad)806: !$OMP& R2,RM2,RM6,RM10,RM12,R1,E1,G1,G2,E1D,G1D,G2D,F_OVER_R,ET) REDUCTION(+:ENERGY,grad)
870: !$OMP DO807: !$OMP DO
871:           do i=1, NC808:           do i=1, NC
872:             C1 = IC(i)809:             C1 = IC(i)
873:             C2 = JC(i)810:             C2 = JC(i)
874:             IF(SBMHESSATOM .EQ. -1 .OR. SBMHESSATOM .EQ. C1 .OR.811:             dx = XSBM(C1) - XSBM(C2)
875:      Q        SBMHESSATOM .EQ. C2)THEN 812:             dy = YSBM(C1) - YSBM(C2)
876:               dx = XSBM(C1) - XSBM(C2)813:             dz = ZSBM(C1) - ZSBM(C2)
877:               dy = YSBM(C1) - YSBM(C2)814:             r2 = dx**2 + dy**2 + dz**2
878:               dz = ZSBM(C1) - ZSBM(C2)815:             J=I*6
879:               r2 = dx**2 + dy**2 + dz**2816:             C13=C1*3
880:               J=I*6817:             C23=C2*3
881:               C13=C1*3818:             if(CONTACTTYPE(I) .eq. 1)then
882:               C23=C2*3819:               ! type 1 is 6-12 interaction
883:               if(CONTACTTYPE(I) .eq. 1)then820:               rm2 = 1.0/r2
884:                 ! type 1 is 6-12 interaction821:               rm6 = rm2**3
885:                 rm2 = 1.0/r2822:               CO1=CONCOEF(J-5)
886:                 rm6 = rm2**3823:               CO2=CONCOEF(J-4)
887:                 CO1=CONCOEF(J-5)824:               ENERGY = ENERGY + rm6*(CO1*rm6-CO2)
888:                 CO2=CONCOEF(J-4)825:               f_over_r = -12.0*rm6*(CO1*rm6-CO2/2.0)*rm2
889:                 ENERGY = ENERGY + rm6*(CO1*rm6-CO2)826:             elseif(CONTACTTYPE(I) .eq. 2)then
890:                 f_over_r = -12.0*rm6*(CO1*rm6-CO2/2.0)*rm2827:               ! type 2 is 10-12 interaction
891:               elseif(CONTACTTYPE(I) .eq. 2)then828:               rm2 = 1.0/r2
892:                 ! type 2 is 10-12 interaction829:               rm10 = rm2**5
893:                 rm2 = 1.0/r2830:               CO1=CONCOEF(J-5)
894:                 rm10 = rm2**5831:               CO2=CONCOEF(J-4)
895:                 CO1=CONCOEF(J-5)832:               ENERGY = ENERGY + rm10*(CO1*rm2-CO2)
896:                 CO2=CONCOEF(J-4)833:               f_over_r = -60.0*rm10*(CO1*rm2/5.0-CO2/6.0)*rm2
897:                 ENERGY = ENERGY + rm10*(CO1*rm2-CO2)834:             elseif(CONTACTTYPE(I) .eq. 5)then
898:                 f_over_r = -60.0*rm10*(CO1*rm2/5.0-CO2/6.0)*rm2835:               ! type 5 is a simple gaussian interaction (no excluded
899:               elseif(CONTACTTYPE(I) .eq. 5)then836:               ! volume)
900:                 ! type 5 is a simple gaussian interaction (no excluded837:               r1 = sqrt(r2)
901:                 ! volume)838:               CO1=CONCOEF(J-5) !A
902:                 r1 = sqrt(r2)839:               CO2=CONCOEF(J-4) !mu
903:                 CO1=CONCOEF(J-5) !A840:               CO3=CONCOEF(J-3) !1/(2*sigma**2)
904:                 CO2=CONCOEF(J-4) !mu841:               E1= -CO1*exp(-CO3*(r1-CO2)**2)
905:                 CO3=CONCOEF(J-3) !1/(2*sigma**2)842:               ENERGY = ENERGY + E1
906:                 E1= -CO1*exp(-CO3*(r1-CO2)**2)843:               f_over_r = -E1*2*CO3*(r1-CO2)/r1
907:                 ENERGY = ENERGY + E1844:             elseif(CONTACTTYPE(I) .eq. 6)then
908:                 f_over_r = -E1*2*CO3*(r1-CO2)/r1845:               ! type 6 is a gaussian interaction with excluded
909:               elseif(CONTACTTYPE(I) .eq. 6)then846:               ! volume included
910:                 ! type 6 is a gaussian interaction with excluded847:               r1 = sqrt(r2)
911:                 ! volume included848:               rm2 = 1.0/r2
912:                 r1 = sqrt(r2)849:               rm12 = rm2**6
913:                 rm2 = 1.0/r2850:               CO1=CONCOEF(J-5) !A
914:                 rm12 = rm2**6851:               CO2=CONCOEF(J-4) !mu
915:                 CO1=CONCOEF(J-5) !A852:               CO3=CONCOEF(J-3) !1/(2*sigma**2)
916:                 CO2=CONCOEF(J-4) !mu853:               CO4=CONCOEF(J-2) !a/A
917:                 CO3=CONCOEF(J-3) !1/(2*sigma**2)854:               E1 =CO1*(1+CO4*rm12)
918:                 CO4=CONCOEF(J-2) !a/A855:               G1 =1-exp(-CO3*(r1-CO2)**2)
919:                 E1 =CO1*(1+CO4*rm12)856:               E1D=-12*CO1*CO4*rm12*rm2
920:                 G1 =1-exp(-CO3*(r1-CO2)**2)857:               G1D=CO3*2*(1-CO2/r1)*exp(-CO3*(r1-CO2)**2)
921:                 E1D=-12*CO1*CO4*rm12*rm2858:               ENERGY = ENERGY + (E1*G1-CO1)
922:                 G1D=CO3*2*(1-CO2/r1)*exp(-CO3*(r1-CO2)**2)859:               f_over_r = E1*G1D + G1*E1D 
923:                 ENERGY = ENERGY + (E1*G1-CO1)860:             elseif(CONTACTTYPE(I) .eq. 7)then
924:                 f_over_r = E1*G1D + G1*E1D 861:               ! type 7 is a double gaussian interaction with excluded
925:               elseif(CONTACTTYPE(I) .eq. 7)then862:               ! volume included
926:                 ! type 7 is a double gaussian interaction with excluded863:               rm2 = 1.0/r2
927:                 ! volume included864:               rm12 = rm2**6
928:                 rm2 = 1.0/r2865:               r1 = sqrt(r2)
929:                 rm12 = rm2**6866:               CO1=CONCOEF(J-5) !A
930:                 r1 = sqrt(r2)867:               CO2=CONCOEF(J-4) !mu1
931:                 CO1=CONCOEF(J-5) !A868:               CO3=CONCOEF(J-3) !1/(2*sigma1**2)
932:                 CO2=CONCOEF(J-4) !mu1869:               CO4=CONCOEF(J-2) !mu2
933:                 CO3=CONCOEF(J-3) !1/(2*sigma1**2)870:               CO5=CONCOEF(J-1) !1/(2*sigma2**2)
934:                 CO4=CONCOEF(J-2) !mu2871:               CO6=CONCOEF(J)   ! a/A (excluded volume amplitude/gaussian amplitude)
935:                 CO5=CONCOEF(J-1) !1/(2*sigma2**2)872:               E1=1+CO6*rm12
936:                 CO6=CONCOEF(J)   ! a/A (excluded volume amplitude/gaussian amplitude)873:               G1=1-exp(-CO3*(r1-CO2)**2)
937:                 E1=1+CO6*rm12874:               G2=1-exp(-CO5*(r1-CO4)**2)
938:                 G1=1-exp(-CO3*(r1-CO2)**2)875:               E1D= -12*CO6*rm12*rm2
939:                 G2=1-exp(-CO5*(r1-CO4)**2)876:               G1D=CO3*2*(1-CO2/r1)*exp(-CO3*(r1-CO2)**2)
940:                 E1D= -12*CO6*rm12*rm2877:               G2D=CO5*2*(1-CO4/r1)*exp(-CO5*(r1-CO4)**2)
941:                 G1D=CO3*2*(1-CO2/r1)*exp(-CO3*(r1-CO2)**2)878:               ENERGY = ENERGY + CO1*( E1 * G1 * G2-1) 
942:                 G2D=CO5*2*(1-CO4/r1)*exp(-CO5*(r1-CO4)**2)879:               f_over_r = CO1*(E1D*G1*G2 
943:                 ENERGY = ENERGY + CO1*( E1 * G1 * G2-1) 880:      Q        +  E1 * G1D * G2
944:                 f_over_r = CO1*(E1D*G1*G2 881:      Q        +  E1 * G1 * G2D)
945:      Q          +  E1 * G1D * G2882:             else
946:      Q          +  E1 * G1 * G2D)883:               WRITE(*,*) 'ERROR: Contact type not recognized:', CONTACTTYPE(I)
947:               else884:               STOP
948:                 WRITE(*,*) 'ERROR: Contact type not recognized:', CONTACTTYPE(I)885:             endif
949:                 STOP886: 
950:               endif887:             grad(C13-2) = grad(C13-2) + f_over_r * dx
951:   888:             grad(C13-1) = grad(C13-1) + f_over_r * dy
952:               grad(C13-2) = grad(C13-2) + f_over_r * dx889:             grad(C13)   = grad(C13)   + f_over_r * dz
953:               grad(C13-1) = grad(C13-1) + f_over_r * dy890:             grad(C23-2) = grad(C23-2) - f_over_r * dx
954:               grad(C13)   = grad(C13)   + f_over_r * dz891:             grad(C23-1) = grad(C23-1) - f_over_r * dy
955:               grad(C23-2) = grad(C23-2) - f_over_r * dx892:             grad(C23)   = grad(C23)   - f_over_r * dz
956:               grad(C23-1) = grad(C23-1) - f_over_r * dy 
957:               grad(C23)   = grad(C23)   - f_over_r * dz 
958:             ENDIF 
959:           enddo893:           enddo
960: !$OMP END DO894: !$OMP END DO
961: !$OMP END PARALLEL895: !$OMP END PARALLEL
962: 896: 
963: 897: 
964:       end898:       end
965: 899: 
966: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^end of SBMContacts^^^^^^^^^^^^^^^^^^^^^^^^^^^900: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^end of SBMContacts^^^^^^^^^^^^^^^^^^^^^^^^^^^
967: 901: 
968: 902: 
970: !* SBMNonContacts computes the forces due to non native contacts      *904: !* SBMNonContacts computes the forces due to non native contacts      *
971: !**********************************************************************905: !**********************************************************************
972: 906: 
973:         subroutine SBMnoncontacts(grad,energy,NATOMS)907:         subroutine SBMnoncontacts(grad,energy,NATOMS)
974:         USE KEY908:         USE KEY
975:         USE SBMDATA909:         USE SBMDATA
976:         implicit NONE910:         implicit NONE
977:         integer I,I3,N,J,AN,NATOMS911:         integer I,I3,N,J,AN,NATOMS
978: 912: 
979:         DOUBLE PRECISION grad(3*NATOMS), energy913:         DOUBLE PRECISION grad(3*NATOMS), energy
980:         DOUBLE PRECISION SAT,SBT,SCT,STTT914:        DOUBLE PRECISION SAT,SBT,SCT,STTT
981: 915: 
982:         integer C1,C2,C13,C23,DINDEX,ii,jj,kk,k,l,iii,jjj,POS916:         integer C1,C2,C13,C23,DINDEX,ii,jj,kk,k,l,iii,jjj,POS
983:         DOUBLE PRECISION  r2, rm2, rm14, f_over_r 917:         DOUBLE PRECISION  r2, rm2, rm14, f_over_r 
984:         INTEGER NTHREADS918:         INTEGER NTHREADS
985: !$    INTEGER  OMP_GET_NUM_THREADS,919: !$    INTEGER  OMP_GET_NUM_THREADS,
986: !$   Q  OMP_GET_THREAD_NUM,OMP_NUM_THREADS920: !$   Q  OMP_GET_THREAD_NUM,OMP_NUM_THREADS
987: 921: 
988:         DOUBLE PRECISION dx,dy,dz922:         DOUBLE PRECISION dx,dy,dz
989:         integer tempN, alpha923:         integer tempN, alpha
990:         double precision Rdiff,Vfunc,Ffunc924:         double precision Rdiff,Vfunc,Ffunc
991:         double precision Rcut2,Rswitch2925:         double precision Rcut2,Rswitch2
992:         ! Ngrid is the number of atoms in that grid point926:         ! Ngrid is the number of atoms in that grid point
993:         ! grid is the array of atoms in each grid927:         ! grid is the array of atoms in each grid
994:         !integer Ngrid,grid928:         integer Ngrid,grid,maxgrid,maxpergrid
995:         integer maxpergrid 
996:         ! number of atoms per grid, max929:         ! number of atoms per grid, max
997:         parameter (maxpergrid=100)930:         parameter (maxpergrid=50)
998:         ! dimensions of grid931:         ! dimensions of grid
999: !        parameter (maxgrid=30)932:         parameter (maxgrid=30)
1000: !        dimension Ngrid(maxgrid,maxgrid,maxgrid),933:         dimension Ngrid(maxgrid,maxgrid,maxgrid),
1001: !     Q  grid(maxgrid,maxgrid,maxgrid,maxpergrid)934:      Q  grid(maxgrid,maxgrid,maxgrid,maxpergrid)
1002:         integer MaxGridX,MaxGridY,MaxGridZ935:         integer MaxGridX,MaxGridY,MaxGridZ
1003:         double precision gridsize,RD1936:         double precision gridsize,RD1
1004:         double precision minX,minY,minZ,maxX,maxY,maxZ937:         double precision minX,minY,minZ,maxX,maxY,maxZ
1005:         integer Xgrid,Ygrid,Zgrid,TID938:         integer Xgrid,Ygrid,Zgrid,TID
1006:         Rdiff=NCcut-NCswitch939:         Rdiff=NCcut-NCswitch
1007:         alpha=12940:         alpha=12
1008:         GRIDSIZE=NCcut*1.01941:         GRIDSIZE=NCcut*1.01
1009:         Rcut2=NCcut**2942:         Rcut2=NCcut**2
1010:         Rswitch2=NCswitch**2943:         Rswitch2=NCswitch**2
1011: 944: 
1012:  
1013:          
1014:         IF(SBMHESSATOM .EQ. -1 .OR. SBMFIRSTDD)THEN  
1015:           SBMFIRSTDD = .FALSE. 
1016:  
1017: ! grid the system945: ! grid the system
1018:         minX=10000000946:         minX=10000000
1019:         minY=10000000947:         minY=10000000
1020:         minZ=10000000948:         minZ=10000000
1021:         949:         
1022:         maxX=-10000000950:         maxX=-10000000
1023:         maxY=-10000000951:         maxY=-10000000
1024:         maxZ=-10000000952:         maxZ=-10000000
1025: 953: 
1026:         do i=1,NATOMS954:         do i=1,NATOMS
1039:                 minZ=ZSBM(i)967:                 minZ=ZSBM(i)
1040:            elseif(ZSBM(i) .gt. maxZ)then968:            elseif(ZSBM(i) .gt. maxZ)then
1041:                 maxZ=ZSBM(i)969:                 maxZ=ZSBM(i)
1042:            endif970:            endif
1043: 971: 
1044:         enddo972:         enddo
1045: 973: 
1046:         maxgridX=int((maxX-minX)/gridsize)+1974:         maxgridX=int((maxX-minX)/gridsize)+1
1047:         maxgridY=int((maxY-minY)/gridsize)+1975:         maxgridY=int((maxY-minY)/gridsize)+1
1048:         maxgridZ=int((maxZ-minZ)/gridsize)+1976:         maxgridZ=int((maxZ-minZ)/gridsize)+1
1049:         if(ALLOCATED(SBMNGRID))THEN977: 
1050:           DEALLOCATE(SBMNGRID)978:         if(maxgridX .ge. maxgrid .or. 
1051:         ENDIF979:      Q  maxgridY .ge. maxgrid .or.
1052:         if(ALLOCATED(SBMGRID))THEN980:      Q  maxgridZ .ge. maxgrid )then
1053:           DEALLOCATE(SBMGRID)981:         write(*,*) 'ERROR: system got too big for grid searching...'
1054:         ENDIF982:         STOP
1055:         ALLOCATE(SBMNGRID(MAXGRIDX,MAXGRIDY,MAXGRIDZ))983:         endif
1056:         ALLOCATE(SBMGRID(MAXGRIDX,MAXGRIDY,MAXGRIDZ,MAXPERGRID))984: 
1057: 985: 
1058: !!  Add a second grid that only includes atoms that are charged986: !!  Add a second grid that only includes atoms that are charged
1059: 987: 
1060:         do i=1,maxgridx988:         do i=1,maxgrid
1061:          do j=1,maxgridy989:          do j=1,maxgrid
1062:           do k=1,maxgridz990:           do k=1,maxgrid
1063:                 SBMNGRID(i,j,k)=0991:                 Ngrid(i,j,k)=0
1064:           enddo992:           enddo
1065:          enddo993:          enddo
1066:         enddo994:         enddo
1067:         do i=1,NATOMS995:         do i=1,NATOMS
1068:                 Xgrid=int((XSBM(i)-minX)/gridsize)+1996:                 Xgrid=int((XSBM(i)-minX)/gridsize)+1
1069:                 Ygrid=int((YSBM(i)-minY)/gridsize)+1997:                 Ygrid=int((YSBM(i)-minY)/gridsize)+1
1070:                 Zgrid=int((ZSBM(i)-minZ)/gridsize)+1998:                 Zgrid=int((ZSBM(i)-minZ)/gridsize)+1
1071:                 SBMNGRID(Xgrid,Ygrid,Zgrid)=SBMNGRID(Xgrid,Ygrid,Zgrid)+1999:                 Ngrid(Xgrid,Ygrid,Zgrid)=Ngrid(Xgrid,Ygrid,Zgrid)+1
1072:                 if(SBMNGRID(Xgrid,Ygrid,Zgrid) .gt. maxpergrid)then1000:                 if(Ngrid(Xgrid,Ygrid,Zgrid) .gt. maxpergrid)then
1073:                         write(*,*) 'ERROR: too many atoms in a grid'1001:                         write(*,*) 'ERROR: too many atoms in a grid'
1074:                         write(*,*) SBMNGRID(Xgrid,Ygrid,Zgrid),Xgrid,Ygrid,Zgrid1002:                         write(*,*) Ngrid(Xgrid,Ygrid,Zgrid),Xgrid,Ygrid,Zgrid
1075:                         STOP1003:                         STOP
1076:                 endif1004:                 endif
1077:                 SBMGRID(Xgrid,Ygrid,Zgrid,SBMNGRID(Xgrid,Ygrid,Zgrid))=i1005:                 grid(Xgrid,Ygrid,Zgrid,Ngrid(Xgrid,Ygrid,Zgrid))=i
1078:         enddo1006:         enddo
1079: 1007: 
1080:         ENDIF 
1081: 1008: 
1082:            tempN=01009:            tempN=0
1083: 1010: 
1084: ! add a second loop that goes over charged atoms1011: ! add a second loop that goes over charged atoms
1085: 1012: 
1086:         IF(SBMHESSATOM .EQ. -1)THEN  
1087:  
1088: !$OMP PARALLEL1013: !$OMP PARALLEL
1089: !$OMP& PRIVATE(i,I3,ii,jj,kk,jjj,C1,C2,C13,C23,r2,rm2,rm14,f_over_r,RD1,dx,dy,dz,1014: !$OMP& PRIVATE(i,I3,ii,jj,kk,jjj,C1,C2,C13,C23,r2,rm2,rm14,f_over_r,RD1,dx,dy,dz,
1090: !$OMP& Xgrid,Ygrid,Zgrid,TID,DINDEX,POS,SAT,SBT,SCT,STTT)  1015: !$OMP& Xgrid,Ygrid,Zgrid,TID,DINDEX,POS,SAT,SBT,SCT,STTT)  
1091: !$OMP& REDUCTION(+:ENERGY,grad)1016: !$OMP& REDUCTION(+:ENERGY,grad)
1092:         TID=01017:         TID=0
1093:         NTHREADS=1;1018:         NTHREADS=1;
1094: !$      TID = OMP_GET_THREAD_NUM()1019: !$      TID = OMP_GET_THREAD_NUM()
1095: !$      NTHREADS = OMP_GET_NUM_THREADS()1020: !$      NTHREADS = OMP_GET_NUM_THREADS()
1096:         do i=1+TID,NATOMS,NTHREADS1021:         do i=1+TID,NATOMS,NTHREADS
1097:  
1098:  
1099:           I3=I*31022:           I3=I*3
1100:           Xgrid=int((XSBM(I)-minX)/gridsize)+11023:           Xgrid=int((XSBM(I)-minX)/gridsize)+1
1101:           Ygrid=int((YSBM(I)-minY)/gridsize)+11024:           Ygrid=int((YSBM(I)-minY)/gridsize)+1
1102:           Zgrid=int((ZSBM(I)-minZ)/gridsize)+11025:           Zgrid=int((ZSBM(I)-minZ)/gridsize)+1
1103:           do ii=XGRID-1,XGRID+11026:           do ii=XGRID-1,XGRID+1
1104:            do jj=YGRID-1,YGRID+11027:            do jj=YGRID-1,YGRID+1
1105:             do kk=ZGRID-1,ZGRID+11028:             do kk=ZGRID-1,ZGRID+1
1106:              if(ii .ge. 1 .and. jj .ge. 1 .and. kk .ge. 1 .and.1029:              if(ii .ge. 1 .and. jj .ge. 1 .and. kk .ge. 1 .and.
1107:      Q            ii .le. MAXGRIDX .and. jj .le. MAXGRIDY .and. kk .le.MAXGRIDZ)then1030:      Q            ii .le. MAXGRIDX .and. jj .le. MAXGRIDY .and. kk .le.MAXGRIDZ)then
1108:               do jjj=1,SBMNGRID(ii,jj,kk)1031:               do jjj=1,Ngrid(ii,jj,kk)
1109:                C2=SBMGRID(ii,jj,kk,jjj)1032:                C2=grid(ii,jj,kk,jjj)
1110:                C23=C2*31033:                C23=C2*3
1111:                DINDEX=C2-I1034:                DINDEX=C2-I
1112:                if(DINDEX .GT. 0)THEN1035:                if(DINDEX .GT. 0)THEN
1113:                 if( DINDEX.GT.MAXSEPSYS.OR.1036:                 if(DINDEX .GT. MAXSEPSYS .OR. (DINDEX .LE. MAXSEPSYS .AND. NNCINC((I-1)*MAXSEP+DINDEX) .eq. 0) )then
1114:      Q             (DINDEX.LE.MAXSEPSYS.AND.NNCINC((I-1)*MAXSEP+DINDEX) .eq. 0))THEN1037:                  if(SBMRBG(I) .NE. SBMRBG(C2))THEN
1115:                    if(SBMRBG(I) .NE. SBMRBG(C2))THEN1038:                   dx = XSBM(I) - XSBM(C2)
1116:                     dx = XSBM(I) - XSBM(C2)1039:                   dy = YSBM(I) - YSBM(C2)
1117:                     dy = YSBM(I) - YSBM(C2)1040:                   dz = ZSBM(I) - ZSBM(C2)
1118:                     dz = ZSBM(I) - ZSBM(C2)1041:                   r2 = dx**2 + dy**2 + dz**2
1119:                     r2 = dx**2 + dy**2 + dz**21042:                   if(r2 .le. Rcut2)then
1120:                     if(r2 .le. Rcut2)then1043:                    POS=(ATOMTYPES(I)-1)*NATOMTYPES+ATOMTYPES(C2)
1121:                      POS=(ATOMTYPES(I)-1)*NATOMTYPES+ATOMTYPES(C2)1044:                    STTT=STT(POS)
1122:                      STTT=STT(POS)1045:                    SCT=SC(POS)
1123:                      SCT=SC(POS)1046:                    rm2 = 1/r2
1124:                      rm2 = 1/r21047:                    rm14 = rm2**7
1125:                      rm14 = rm2**71048:                    energy=energy+STTT*rm2**6+SCT
1126:                      energy=energy+STTT*rm2**6+SCT1049:                    f_over_r=-STTT*12.0*rm14
1127:                      f_over_r=-STTT*12.0*rm141050:                    RD1=sqrt(r2)-NCswitch
1128:                      RD1=sqrt(r2)-NCswitch1051:                    if(r2 .gt. Rswitch2)then
1129:                      if(r2 .gt. Rswitch2)then1052:                     SAT=SA(POS)
1130:                       SAT=SA(POS)1053:                     SBT=SB(POS)
1131:                       SBT=SB(POS)1054:                     f_over_r=f_over_r+(SAT*RD1**2+SBT*RD1**3)*sqrt(rm2)
1132:                       f_over_r=f_over_r+(SAT*RD1**2+SBT*RD1**3)*sqrt(rm2)1055:                     energy=energy+SAT/3.0*RD1**3+SBT/4.0*RD1**4
1133:                       energy=energy+SAT/3.0*RD1**3+SBT/4.0*RD1**41056:                    elseif(r2 .lt. Rswitch2)then
1134:                      endif1057:                   ! normal repulsive term
1135:                      grad(I3-2) = grad(I3-2) + f_over_r * dx1058:                    else
1136:                      grad(I3-1) = grad(I3-1) + f_over_r * dy1059:                   ! things should have fallen in one of the previous two...
1137:                      grad(I3)   = grad(I3)   + f_over_r * dz1060:                     write(*,*) 'something went wrong with switching function'
1138:                      grad(C23-2) = grad(C23-2) - f_over_r * dx1061:                    endif
1139:                      grad(C23-1) = grad(C23-1) - f_over_r * dy1062:                    grad(I3-2) = grad(I3-2) + f_over_r * dx
1140:                      grad(C23)   = grad(C23)   - f_over_r * dz1063:                    grad(I3-1) = grad(I3-1) + f_over_r * dy
 1064:                    grad(I3)   = grad(I3)   + f_over_r * dz
 1065:                    grad(C23-2) = grad(C23-2) - f_over_r * dx
 1066:                    grad(C23-1) = grad(C23-1) - f_over_r * dy
 1067:                    grad(C23)   = grad(C23)   - f_over_r * dz
1141:                   ENDIF1068:                   ENDIF
1142:                   endif1069:                   endif
1143:                 endif1070:                 endif
1144:                endif1071:                endif
1145:               enddo1072:               enddo
1146:              endif1073:              endif
1147:             enddo1074:             enddo
1148:            enddo1075:            enddo
1149:           enddo1076:           enddo
1150:          enddo1077:          enddo
1151: 1078: 
1152: !$OMP DO1079: !$OMP DO
1153:        DO I=1,NEXCLUSIONS1080:        DO I=1,NEXCLUSIONS
1154:          C1=NNEXL1(I) 1081:          C1=NNEXL1(I) 
1155:          C2=NNEXL2(I)1082:          C2=NNEXL2(I)
1156:            C13=C1*31083:          C13=C1*3
1157:            C23=C2*31084:          C23=C2*3
1158:            dx = XSBM(C1) - XSBM(C2)1085:          dx = XSBM(C1) - XSBM(C2)
1159:            dy = YSBM(C1) - YSBM(C2)1086:          dy = YSBM(C1) - YSBM(C2)
1160:            dz = ZSBM(C1) - ZSBM(C2)1087:          dz = ZSBM(C1) - ZSBM(C2)
1161:            r2 = dx**2 + dy**2 + dz**21088:          r2 = dx**2 + dy**2 + dz**2
1162:            if(r2 .le. Rcut2)then1089:          if(r2 .le. Rcut2)then
1163:              POS=(ATOMTYPES(C1)-1)*NATOMTYPES+ATOMTYPES(C2)1090:            POS=(ATOMTYPES(C1)-1)*NATOMTYPES+ATOMTYPES(C2)
1164:              STTT=STT(POS)1091:            STTT=STT(POS)
1165:              SCT=SC(POS)1092:            SCT=SC(POS)
1166:    1093:  
1167:              rm2 = 1/r21094:            rm2 = 1/r2
1168:              rm14 = rm2**71095:            rm14 = rm2**7
1169:               energy=energy-STTT*rm2**6-SCT1096:             energy=energy-STTT*rm2**6-SCT
1170:               f_over_r=-STTT*12.0*rm141097:             f_over_r=-STTT*12.0*rm14
1171:               RD1=sqrt(r2)-NCswitch1098:             RD1=sqrt(r2)-NCswitch
1172:             if(r2 .gt. Rswitch2)then1099:           if(r2 .gt. Rswitch2)then
1173:              SAT=SA(POS)1100:            SAT=SA(POS)
1174:              SBT=SB(POS)1101:            SBT=SB(POS)
1175:              f_over_r=f_over_r+(SAT*RD1**2+SBT*RD1**3)*sqrt(rm2)1102:            f_over_r=f_over_r+(SAT*RD1**2+SBT*RD1**3)*sqrt(rm2)
1176:              energy=energy-SAT/3.0*RD1**3-SBT/4.0*RD1**41103:            energy=energy-SAT/3.0*RD1**3-SBT/4.0*RD1**4
1177:             endif1104: !          elseif(r2 .lt. Rswitch2)then
1178:             grad(C13-2) = grad(C13-2) - f_over_r * dx1105:          ! normal repulsive term
1179:             grad(C13-1) = grad(C13-1) - f_over_r * dy1106: !          else
1180:             grad(C13)   = grad(C13)   - f_over_r * dz1107:          ! things should have fallen in one of the previous two...
1181:             grad(C23-2) = grad(C23-2) + f_over_r * dx1108: !           write(*,*) 'something went wrong with switching function'
1182:             grad(C23-1) = grad(C23-1) + f_over_r * dy1109:           endif
1183:             grad(C23)   = grad(C23)   + f_over_r * dz1110:           grad(C13-2) = grad(C13-2) - f_over_r * dx
1184:            endif1111:           grad(C13-1) = grad(C13-1) - f_over_r * dy
 1112:           grad(C13)   = grad(C13)   - f_over_r * dz
 1113:           grad(C23-2) = grad(C23-2) + f_over_r * dx
 1114:           grad(C23-1) = grad(C23-1) + f_over_r * dy
 1115:           grad(C23)   = grad(C23)   + f_over_r * dz
 1116:          endif
1185:       ENDDO1117:       ENDDO
1186: !$OMP END DO1118: !$OMP END DO
1187: !$OMP END PARALLEL1119: !$OMP END PARALLEL
1188:   
1189:  
1190:          ELSE  ! SBMHESSATOM IS DEFINED 
1191:  
1192:           I=SBMHESSATOM 
1193:  
1194:           I3=I*3 
1195:           Xgrid=int((XSBM(I)-minX)/gridsize)+1 
1196:           Ygrid=int((YSBM(I)-minY)/gridsize)+1 
1197:           Zgrid=int((ZSBM(I)-minZ)/gridsize)+1 
1198:           do ii=XGRID-1,XGRID+1 
1199:            do jj=YGRID-1,YGRID+1 
1200:             do kk=ZGRID-1,ZGRID+1 
1201:              if(ii .ge. 1 .and. jj .ge. 1 .and. kk .ge. 1 .and. 
1202:      Q            ii .le. MAXGRIDX .and. jj .le. MAXGRIDY .and. kk .le.MAXGRIDZ)then 
1203:               do jjj=1,SBMNGRID(ii,jj,kk) 
1204:                C2=SBMGRID(ii,jj,kk,jjj) 
1205:                DINDEX=C2-I 
1206: !               if(DINDEX .GT. 0 .OR. SBMHESSATOM .EQ. I)THEN 
1207:                 if( C2 .NE. I)THEN 
1208:                  C23=C2*3 
1209:                 if( DINDEX.GT.MAXSEPSYS.OR. 
1210:      Q             (DINDEX.LE.MAXSEPSYS.AND.NNCINC((I-1)*MAXSEP+DINDEX) .eq. 0) .OR. 
1211:      Q             (DINDEX.LT.-MAXSEPSYS) .OR. 
1212:      Q             (DINDEX.GE.-MAXSEPSYS.AND.NNCINC((C2-1)*MAXSEP-DINDEX) .eq. 0)  )then 
1213:                    if(SBMRBG(I) .NE. SBMRBG(C2))THEN 
1214:                     dx = XSBM(I) - XSBM(C2) 
1215:                     dy = YSBM(I) - YSBM(C2) 
1216:                     dz = ZSBM(I) - ZSBM(C2) 
1217:                     r2 = dx**2 + dy**2 + dz**2 
1218:                     if(r2 .le. Rcut2)then 
1219:                      POS=(ATOMTYPES(I)-1)*NATOMTYPES+ATOMTYPES(C2) 
1220:                      STTT=STT(POS) 
1221:                      SCT=SC(POS) 
1222:                      rm2 = 1/r2 
1223:                      rm14 = rm2**7 
1224:                      energy=energy+STTT*rm2**6+SCT 
1225:                      f_over_r=-STTT*12.0*rm14 
1226:                      RD1=sqrt(r2)-NCswitch 
1227:                      if(r2 .gt. Rswitch2)then 
1228:                       SAT=SA(POS) 
1229:                       SBT=SB(POS) 
1230:                       f_over_r=f_over_r+(SAT*RD1**2+SBT*RD1**3)*sqrt(rm2) 
1231:                       energy=energy+SAT/3.0*RD1**3+SBT/4.0*RD1**4 
1232:                      elseif(r2 .lt. Rswitch2)then 
1233:                     ! normal repulsive term 
1234:                      else 
1235:                     ! things should have fallen in one of the previous two... 
1236:                       write(*,*) 'something went wrong with switching function' 
1237:                      endif 
1238:                      grad(I3-2) = grad(I3-2) + f_over_r * dx 
1239:                      grad(I3-1) = grad(I3-1) + f_over_r * dy 
1240:                      grad(I3)   = grad(I3)   + f_over_r * dz 
1241:                      grad(C23-2) = grad(C23-2) - f_over_r * dx 
1242:                      grad(C23-1) = grad(C23-1) - f_over_r * dy 
1243:                      grad(C23)   = grad(C23)   - f_over_r * dz 
1244:                   ENDIF 
1245:                   endif 
1246:                 endif 
1247:                endif 
1248:               enddo 
1249:              endif 
1250:             enddo 
1251:            enddo 
1252:           enddo 
1253: !!         enddo 
1254:  
1255:  
1256:        DO I=1,NEXCLUSIONS 
1257:          C1=NNEXL1(I)  
1258:          C2=NNEXL2(I) 
1259:          IF(SBMHESSATOM .EQ. C1 .OR. SBMHESSATOM .EQ. C2)THEN  
1260:            C13=C1*3 
1261:            C23=C2*3 
1262:            dx = XSBM(C1) - XSBM(C2) 
1263:            dy = YSBM(C1) - YSBM(C2) 
1264:            dz = ZSBM(C1) - ZSBM(C2) 
1265:            r2 = dx**2 + dy**2 + dz**2 
1266:            if(r2 .le. Rcut2)then 
1267:              POS=(ATOMTYPES(C1)-1)*NATOMTYPES+ATOMTYPES(C2) 
1268:              STTT=STT(POS) 
1269:              SCT=SC(POS) 
1270:     
1271:              rm2 = 1/r2 
1272:              rm14 = rm2**7 
1273:               energy=energy-STTT*rm2**6-SCT 
1274:               f_over_r=-STTT*12.0*rm14 
1275:               RD1=sqrt(r2)-NCswitch 
1276:             if(r2 .gt. Rswitch2)then 
1277:              SAT=SA(POS) 
1278:              SBT=SB(POS) 
1279:              f_over_r=f_over_r+(SAT*RD1**2+SBT*RD1**3)*sqrt(rm2) 
1280:              energy=energy-SAT/3.0*RD1**3-SBT/4.0*RD1**4 
1281:             endif 
1282:             grad(C13-2) = grad(C13-2) - f_over_r * dx 
1283:             grad(C13-1) = grad(C13-1) - f_over_r * dy 
1284:             grad(C13)   = grad(C13)   - f_over_r * dz 
1285:             grad(C23-2) = grad(C23-2) + f_over_r * dx 
1286:             grad(C23-1) = grad(C23-1) + f_over_r * dy 
1287:             grad(C23)   = grad(C23)   + f_over_r * dz 
1288:            endif 
1289:          ENDIF 
1290:        ENDDO 
1291:       ENDIF 
1292:  
1293:  
1294:       END1120:       END
1295: 1121: 
1296: !^^^^^^^^^^^^^^^^^^^^^^^^^^^End of SBMNonContacts^^^^^^^^^^^^^^^^^^^^^^^^^1122: !^^^^^^^^^^^^^^^^^^^^^^^^^^^End of SBMNonContacts^^^^^^^^^^^^^^^^^^^^^^^^^
1297: 1123: 
1298: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<1124: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1299: ! FUNCTIONS for DH Calculations1125: ! FUNCTIONS for DH Calculations
1300: !*****************************************************1126: !*****************************************************
1301: 1127: 
1302:       DOUBLE PRECISION FUNCTION SBMDHV(kappa,r)1128:       DOUBLE PRECISION FUNCTION SBMDHV(kappa,r)
1303:       USE KEY1129:       USE KEY
1333:         INTEGER NTHREADS,  OMP_GET_NUM_THREADS,1159:         INTEGER NTHREADS,  OMP_GET_NUM_THREADS,
1334:      Q  OMP_GET_THREAD_NUM,OMP_NUM_THREADS1160:      Q  OMP_GET_THREAD_NUM,OMP_NUM_THREADS
1335:         DOUBLE PRECISION dx,dy,dz1161:         DOUBLE PRECISION dx,dy,dz
1336:         integer tempN, alpha1162:         integer tempN, alpha
1337:         double precision diff,Vfunc,Ffunc1163:         double precision diff,Vfunc,Ffunc
1338:         double precision DHswitch2,DHcut21164:         double precision DHswitch2,DHcut2
1339:         ! Ngrid is the number of atoms in that grid point1165:         ! Ngrid is the number of atoms in that grid point
1340:         ! grid is the array of atoms in each grid1166:         ! grid is the array of atoms in each grid
1341:         integer Ngrid,grid,maxgrid,maxpergrid1167:         integer Ngrid,grid,maxgrid,maxpergrid
1342:         ! number of atoms per grid, max1168:         ! number of atoms per grid, max
1343:         parameter (maxpergrid=100)1169:         parameter (maxpergrid=50)
1344:         ! dimensions of grid1170:         ! dimensions of grid
1345:         parameter (maxgrid=15)1171:         parameter (maxgrid=15)
1346:         dimension Ngrid(maxgrid,maxgrid,maxgrid),1172:         dimension Ngrid(maxgrid,maxgrid,maxgrid),
1347:      Q  grid(maxgrid,maxgrid,maxgrid,maxpergrid)1173:      Q  grid(maxgrid,maxgrid,maxgrid,maxpergrid)
1348:         integer MaxGridX,MaxGridY,MaxGridZ1174:         integer MaxGridX,MaxGridY,MaxGridZ
1349:         double precision gridsize,RD11175:         double precision gridsize,RD1
1350:         double precision minX,minY,minZ,maxX,maxY,maxZ1176:         double precision minX,minY,minZ,maxX,maxY,maxZ
1351:         integer Xgrid,Ygrid,Zgrid,TID1177:         integer Xgrid,Ygrid,Zgrid,TID
1352:         double precision C2T1178:         double precision C2T
1353: 1179: 
1450:           Ygrid=int((YSBM(C1)-minY)/gridsize)+11276:           Ygrid=int((YSBM(C1)-minY)/gridsize)+1
1451:           Zgrid=int((ZSBM(C1)-minZ)/gridsize)+11277:           Zgrid=int((ZSBM(C1)-minZ)/gridsize)+1
1452:         1278:         
1453:           do ii=XGRID-1,XGRID+11279:           do ii=XGRID-1,XGRID+1
1454:            do jj=YGRID-1,YGRID+11280:            do jj=YGRID-1,YGRID+1
1455:             do kk=ZGRID-1,ZGRID+11281:             do kk=ZGRID-1,ZGRID+1
1456:              if(ii .ge. 1 .and. jj .ge. 1 .and. kk .ge. 1 .and.1282:              if(ii .ge. 1 .and. jj .ge. 1 .and. kk .ge. 1 .and.
1457:      Q            ii .le. MAXGRIDX .and. jj .le. MAXGRIDY .and. kk .le.MAXGRIDZ)then1283:      Q            ii .le. MAXGRIDX .and. jj .le. MAXGRIDY .and. kk .le.MAXGRIDZ)then
1458:               do jjj=1,Ngrid(ii,jj,kk)1284:               do jjj=1,Ngrid(ii,jj,kk)
1459:                C2=grid(ii,jj,kk,jjj)1285:                C2=grid(ii,jj,kk,jjj)
1460:                IF(SBMHESSATOM .EQ. -1 .OR. SBMHESSATOM .EQ. C1 .OR.1286:                C23=C2*3
1461:      Q          SBMHESSATOM .EQ. C2)THEN 1287:                if(C2-C1 .gt. 0)then
1462:                 C23=C2*31288:                 dx = XSBM(C1) - XSBM(C2)
1463:                 if(C2-C1 .gt. 0)then1289:                 dy = YSBM(C1) - YSBM(C2)
1464:                  dx = XSBM(C1) - XSBM(C2)1290:                 dz = ZSBM(C1) - ZSBM(C2)
1465:                  dy = YSBM(C1) - YSBM(C2)1291:                 r2 = dx**2 + dy**2 + dz**2
1466:                  dz = ZSBM(C1) - ZSBM(C2)1292:                 if(r2 .le. DHcut2)then
1467:                  r2 = dx**2 + dy**2 + dz**21293:                  C2T=SBMCHARGE(C1)*SBMCHARGE(C2)
1468:                    if(r2 .le. DHcut2)then1294:         ! add force evaluation now
1469:                   C2T=SBMCHARGE(C1)*SBMCHARGE(C2)1295:                  rm2 = 1/r2
1470:           ! add force evaluation now1296:                  r1=sqrt(r2)
1471:                   rm2 = 1/r21297:                  energy=energy+PREFACTOR*C2T*SBMDHV(kappa,r1)
1472:                   r1=sqrt(r2)1298:                  f_over_r=PREFACTOR*C2T*SBMDHVP(kappa,r1)
1473:                     energy=energy+PREFACTOR*C2T*SBMDHV(kappa,r1)1299:                  if(r2 .gt. DHswitch2)then
1474:                     f_over_r=PREFACTOR*C2T*SBMDHVP(kappa,r1)1300:                   RD1=r1-DHswitch
1475:                       if(r2 .gt. DHswitch2)then1301:                   f_over_r=(f_over_r+C2T*(2*COEF2*RD1+3*COEF3*RD1**2))
1476:                        RD1=r1-DHswitch1302:                   energy=energy+COEF2*RD1**2+COEF3*RD1**3
1477:                        f_over_r=(f_over_r+C2T*(2*COEF2*RD1+3*COEF3*RD1**2))1303:                  elseif(r2 .lt. DHswitch2)then
1478:                        energy=energy+COEF2*RD1**2+COEF3*RD1**31304:                 ! normal DH term
1479:                       elseif(r2 .lt. DHswitch2)then1305:                  else
1480:                    ! normal DH term1306:                 ! things should have fallen in one of the previous two...
1481:                       else1307:                   write(*,*) 'something went wrong with DH switching function'
1482:                    ! things should have fallen in one of the previous two...1308:                   STOP
1483:                        write(*,*) 'something went wrong with DH switching function'1309:                  endif
1484:                     STOP1310:                  f_over_r=f_over_r*sqrt(rm2)
1485:                       endif1311:                  grad(C13-2) = grad(C13-2) + f_over_r * dx
1486:                    f_over_r=f_over_r*sqrt(rm2)1312:                  grad(C13-1) = grad(C13-1) + f_over_r * dy
1487:                    grad(C13-2) = grad(C13-2) + f_over_r * dx1313:                  grad(C13)   = grad(C13)   + f_over_r * dz
1488:                    grad(C13-1) = grad(C13-1) + f_over_r * dy1314:                  grad(C23-2) = grad(C23-2) - f_over_r * dx
1489:                    grad(C13)   = grad(C13)   + f_over_r * dz1315:                  grad(C23-1) = grad(C23-1) - f_over_r * dy
1490:                    grad(C23-2) = grad(C23-2) - f_over_r * dx1316:                  grad(C23)   = grad(C23)   - f_over_r * dz
1491:                    grad(C23-1) = grad(C23-1) - f_over_r * dy1317: 
1492:                    grad(C23)   = grad(C23)   - f_over_r * dz 
1493:                   ENDIF 
1494:                   endif1318:                   endif
1495:                 endif1319:                 endif
1496:                enddo1320:                enddo
1497:                   endif1321:                   endif
1498:              enddo1322:              enddo
1499:             enddo1323:             enddo
1500:            enddo1324:            enddo
1501:           enddo1325:           enddo
1502: 1326: 
1503: !!! Add routine to subtract exclusion interactions1327: !!! Add routine to subtract exclusion interactions
1504: !$OMP DO1328: !$OMP DO
1505:        DO I=1,NEXCLUSIONS1329:        DO I=1,NEXCLUSIONS
1506:          C1=NNEXL1(I) 1330:          C1=NNEXL1(I) 
1507:          C2=NNEXL2(I)1331:          C2=NNEXL2(I)
1508:          IF(SBMHESSATOM .EQ. -1 .OR. SBMHESSATOM .EQ. C1 .OR.1332:          C13=C1*3
1509:      Q     SBMHESSATOM .EQ. C2)THEN 1333:          C23=C2*3
1510:            C13=C1*31334:          dx = XSBM(C1) - XSBM(C2)
1511:            C23=C2*31335:          dy = YSBM(C1) - YSBM(C2)
1512:            dx = XSBM(C1) - XSBM(C2)1336:          dz = ZSBM(C1) - ZSBM(C2)
1513:            dy = YSBM(C1) - YSBM(C2)1337:          r2 = dx**2 + dy**2 + dz**2
1514:            dz = ZSBM(C1) - ZSBM(C2)1338:          if(r2 .le. DHcut2)then
1515:            r2 = dx**2 + dy**2 + dz**21339:           C2T=SBMCHARGE(C1)*SBMCHARGE(C2)
1516:            if(r2 .le. DHcut2)then1340:           rm2 = 1/r2
1517:             C2T=SBMCHARGE(C1)*SBMCHARGE(C2)1341:           r1=sqrt(r2)
1518:             rm2 = 1/r21342:           energy=energy-PREFACTOR*C2T*SBMDHV(kappa,r1)
1519:             r1=sqrt(r2)1343:           f_over_r=-PREFACTOR*C2T*SBMDHVP(kappa,r1)
1520:             energy=energy-PREFACTOR*C2T*SBMDHV(kappa,r1)1344:           if(r2 .gt. DHswitch2)then
1521:             f_over_r=-PREFACTOR*C2T*SBMDHVP(kappa,r1)1345:            RD1=r1-DHswitch
1522:             if(r2 .gt. DHswitch2)then1346:            f_over_r=f_over_r-C2T*(2*COEF2*RD1+3*COEF3*RD1**2)
1523:              RD1=r1-DHswitch1347:            energy=energy-COEF2*RD1**2+COEF3*RD1**3
1524:              f_over_r=f_over_r-C2T*(2*COEF2*RD1+3*COEF3*RD1**2)1348:           elseif(r2 .lt. DHswitch2)then
1525:              energy=energy-COEF2*RD1**2+COEF3*RD1**31349:          ! normal DH term
1526:             elseif(r2 .lt. DHswitch2)then1350:           else
1527:            ! normal DH term1351:          ! things should have fallen in one of the previous two...
1528:             else1352:            write(*,*) 'something went wrong with DH switching function'
1529:            ! things should have fallen in one of the previous two...1353:           endif
1530:             write(*,*) 'something went wrong with DH switching function'1354:           f_over_r=f_over_r*sqrt(rm2)
1531:             endif1355:           grad(C13-2) = grad(C13-2) + f_over_r * dx
1532:             f_over_r=f_over_r*sqrt(rm2)1356:           grad(C13-1) = grad(C13-1) + f_over_r * dy
1533:             grad(C13-2) = grad(C13-2) + f_over_r * dx1357:           grad(C13)   = grad(C13)   + f_over_r * dz
1534:             grad(C13-1) = grad(C13-1) + f_over_r * dy1358:           grad(C23-2) = grad(C23-2) - f_over_r * dx
1535:             grad(C13)   = grad(C13)   + f_over_r * dz1359:           grad(C23-1) = grad(C23-1) - f_over_r * dy
1536:             grad(C23-2) = grad(C23-2) - f_over_r * dx1360:           grad(C23)   = grad(C23)   - f_over_r * dz
1537:             grad(C23-1) = grad(C23-1) - f_over_r * dy 
1538:             grad(C23)   = grad(C23)   - f_over_r * dz 
1539:            ENDIF 
1540:          ENDIF1361:          ENDIF
1541:        ENDDO1362:        ENDDO
1542: !$OMP END DO1363: !$OMP END DO
1543: !$OMP END PARALLEL1364: !$OMP END PARALLEL
1544: 1365: 
1545:       END1366:       END
1546: 1367: 
1547: !^^^^^^^^^^^^^^^^^^^^^^^^^^^End of SBMDHELEC^^^^^^^^^^^^^^^^^^^^^^^^^1368: !^^^^^^^^^^^^^^^^^^^^^^^^^^^End of SBMDHELEC^^^^^^^^^^^^^^^^^^^^^^^^^
1548: 1369: 
1549: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<1370: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1559:       DOUBLE PRECISION grad(3*NATOMS), energy,ET1380:       DOUBLE PRECISION grad(3*NATOMS), energy,ET
1560:       DOUBLE PRECISION DX,DY,DZ,K1,K2,K3,K4,K5,K61381:       DOUBLE PRECISION DX,DY,DZ,K1,K2,K3,K4,K5,K6
1561: 1382: 
1562: 1383: 
1563:       if(SBMPRN .eq. 0) RETURN1384:       if(SBMPRN .eq. 0) RETURN
1564: 1385: 
1565: !$OMP PARALLEL PRIVATE(I,J,J3,DX,DY,DZ,K1,K2,K3,K4,K5,K6)REDUCTION(+:energy,grad)1386: !$OMP PARALLEL PRIVATE(I,J,J3,DX,DY,DZ,K1,K2,K3,K4,K5,K6)REDUCTION(+:energy,grad)
1566: !$OMP DO1387: !$OMP DO
1567:               do I=1,SBMPRN1388:               do I=1,SBMPRN
1568:            J=SBMPRI(I)1389:            J=SBMPRI(I)
1569:  
1570:          IF(SBMHESSATOM .EQ. -1 .OR. SBMHESSATOM .EQ. J )THEN  
1571:            I1=(I-1)*61390:            I1=(I-1)*6
1572:            I2=(I-1)*31391:            I2=(I-1)*3
1573:            J3=J*31392:            J3=J*3
1574:            DX=XSBM(J)-SBMPRX(I2+1)1393:            DX=XSBM(J)-SBMPRX(I2+1)
1575:            DY=YSBM(J)-SBMPRX(I2+2)1394:            DY=YSBM(J)-SBMPRX(I2+2)
1576:            DZ=ZSBM(J)-SBMPRX(I2+3)1395:            DZ=ZSBM(J)-SBMPRX(I2+3)
1577:            K1=SBMPRK(I1+1)1396:            K1=SBMPRK(I1+1)
1578:            K2=SBMPRK(I1+2)1397:            K2=SBMPRK(I1+2)
1579:            K3=SBMPRK(I1+3)1398:            K3=SBMPRK(I1+3)
1580:            K4=SBMPRK(I1+4)1399:            K4=SBMPRK(I1+4)
1581:            K5=SBMPRK(I1+5)1400:            K5=SBMPRK(I1+5)
1582:            K6=SBMPRK(I1+6)1401:            K6=SBMPRK(I1+6)
1583:            1402:            
1584:            ENERGY = ENERGY +0.5*(K1*DX**2+K2*DY**2+K3*DZ**2) +1403:            ENERGY = ENERGY +0.5*(K1*DX**2+K2*DY**2+K3*DZ**2) +
1585:      Q     K4*DX*DY+K5*DX*DZ+K6*DY*DZ1404:      Q     K4*DX*DY+K5*DX*DZ+K6*DY*DZ
1586: 1405: 
1587:            grad(J3-2) = grad(J3-2) + K1*DX + K4*DY + K5*DZ 1406:            grad(J3-2) = grad(J3-2) + K1*DX + K4*DY + K5*DZ 
1588:            grad(J3-1) = grad(J3-1) + K2*DY + K4*DX + K6*DZ1407:            grad(J3-1) = grad(J3-1) + K2*DY + K4*DX + K6*DZ
1589:            grad(J3)   = grad(J3)   + K3*DZ + K5*DX + K6*DY1408:            grad(J3)   = grad(J3)   + K3*DZ + K5*DX + K6*DY
1590:          ENDIF 
1591:         enddo1409:         enddo
1592: !$OMP END DO1410: !$OMP END DO
1593: !$OMP END PARALLEL1411: !$OMP END PARALLEL
1594: 1412: 
1595:       END1413:       END
1596: 1414: 
1597: !^^^^^^^^^^^^^^^^^^^^^^^^^^^End of SBMPR^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^1415: !^^^^^^^^^^^^^^^^^^^^^^^^^^^End of SBMPR^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1598: 1416: 


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0