hdiff output

r30664/SBM.f 2016-07-06 15:35:47.475385755 +0100 r30663/SBM.f 2016-07-06 15:35:47.819390407 +0100
334: 334: 
335: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<335: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
336: !* SBMBonds  computes the hookean force and energy between chosen atoms *336: !* SBMBonds  computes the hookean force and energy between chosen atoms *
337: !***********************************************************************337: !***********************************************************************
338: 338: 
339:       subroutine SBMBonds(x,y,z,grad,energy, natoms,Ib1, Ib2,Rb, bK,NBA)339:       subroutine SBMBonds(x,y,z,grad,energy, natoms,Ib1, Ib2,Rb, bK,NBA)
340:       USE KEY340:       USE KEY
341:       implicit NONE341:       implicit NONE
342:       integer I2,J2,outE,I,N,J,NATOMS,NBA,C1342:       integer I2,J2,outE,I,N,J,NATOMS,NBA,C1
343:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS), grad(3*NATOMS),343:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS), grad(3*NATOMS),
344:      Q energy344:      Q energy,ETEMP
345:       DOUBLE PRECISION r2, f, r1345:       DOUBLE PRECISION r2, f, r1
346:       DOUBLE PRECISION dx,dy,dz346:       DOUBLE PRECISION dx,dy,dz
347:       DOUBLE PRECISION Rb(NBA), bK(NBA)347:       DOUBLE PRECISION Rb(NBA), bK(NBA),GRADOMP(3*NATOMS)
348:       INTEGER Ib1(NBA), Ib2(NBA)348:       INTEGER Ib1(NBA), Ib2(NBA)
349: 349: 
350: !$OMP PARALLEL PRIVATE(I,I2,J2,DX,DY,DZ,R1,R2,F)REDUCTION(+:ENERGY,grad)350:       do I=1,NATOMS
 351:          gradOMP(I*3)=0
 352:          gradOMP(I*3-1)=0
 353:          gradOMP(I*3-2)=0
 354:       enddo
 355:         etemp=0
 356: !$OMP PARALLEL PRIVATE(I,I2,J2,DX,DY,DZ,R1,R2,F)REDUCTION(+:etemp,gradOMP)
351: !$OMP DO357: !$OMP DO
352:         DO I=1, nBA358:         DO I=1, nBA
353:            I2 = Ib1(I)359:            I2 = Ib1(I)
354:            J2 = Ib2(I)360:            J2 = Ib2(I)
355:            dx = X(I2) - X(J2)361:            dx = X(I2) - X(J2)
356:            dy = Y(I2) - Y(J2)362:            dy = Y(I2) - Y(J2)
357:            dz = Z(I2) - Z(J2)363:            dz = Z(I2) - Z(J2)
358:            r2 = dx**2 + dy**2 + dz**2364:            r2 = dx**2 + dy**2 + dz**2
359:            r1 = sqrt(r2)365:            r1 = sqrt(r2)
360:            ENERGY = ENERGY + bk(I)*(r1-Rb(I))**2/2.0366:            etemp = etemp + bk(I)*(r1-Rb(I))**2/2.0
361:            f = -bk(I)*(r1-Rb(I))/r1367:            f = -bk(I)*(r1-Rb(I))/r1
362:            grad(I2*3-2) = grad(I2*3-2) - f * dx368:            gradOMP(I2*3-2) = gradOMP(I2*3-2) - f * dx
363:            grad(I2*3-1) = grad(I2*3-1) - f * dy369:            gradOMP(I2*3-1) = gradOMP(I2*3-1) - f * dy
364:            grad(I2*3)   = grad(I2*3)   - f * dz370:            gradOMP(I2*3)   = gradOMP(I2*3)   - f * dz
365:            grad(J2*3-2) = grad(J2*3-2) + f * dx371:            gradOMP(J2*3-2) = gradOMP(J2*3-2) + f * dx
366:            grad(J2*3-1) = grad(J2*3-1) + f * dy372:            gradOMP(J2*3-1) = gradOMP(J2*3-1) + f * dy
367:            grad(J2*3)   = grad(J2*3)   + f * dz373:            gradOMP(J2*3)   = gradOMP(J2*3)   + f * dz
368:        ENDDO 374:        ENDDO 
369: !$OMP END DO375: !$OMP END DO
370: !$OMP END PARALLEL376: !$OMP END PARALLEL
 377:         energy=energy+etemp 
 378:         do C1=1,NATOMS
 379:               grad(C1*3-2) = grad(C1*3-2)+ gradOMP(C1*3-2) 
 380:               grad(C1*3-1) = grad(C1*3-1)+ gradOMP(C1*3-1) 
 381:               grad(C1*3)   = grad(C1*3)  + gradOMP(C1*3)   
 382:         enddo
371: 383: 
372: 384: 
373:       END385:       END
374: 386: 
375: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^END OF SBMBONDS^^^^^^^^^^^^^^^^^^^^^^^^^^^^^387: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^END OF SBMBONDS^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
376: 388: 
377: 389: 
378: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<390: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
379: !* SBMANGL  computes the Force due to the bond angles                   *391: !* SBMANGL  computes the Force due to the bond angles                   *
380: !* This code is modeled after how AMBER performs angle forces           *392: !* This code is modeled after how AMBER performs angle forces           *
384:      Q ANTC, Tk,NTA)396:      Q ANTC, Tk,NTA)
385:       USE KEY397:       USE KEY
386:       implicit NONE398:       implicit NONE
387:       integer NATOMS,C1399:       integer NATOMS,C1
388:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS), grad(3*NATOMS),400:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS), grad(3*NATOMS),
389:      Q energy401:      Q energy
390:       integer NTA402:       integer NTA
391:       LOGICAL SKIP,NOCRST403:       LOGICAL SKIP,NOCRST
392: 404: 
393:       DOUBLE PRECISION RIJ,RKJ,RIK,ANT,XIJ,YIJ,405:       DOUBLE PRECISION RIJ,RKJ,RIK,ANT,XIJ,YIJ,
394:      + ZIJ,XKJ,YKJ,ZKJ, DF406:      + ZIJ,XKJ,YKJ,ZKJ, DF,etemp,gradOMP
 407:       DIMENSION gradOMP(3*NATOMS)
395:       DOUBLE PRECISION CT0, CT1, CT2, DA, ST, 408:       DOUBLE PRECISION CT0, CT1, CT2, DA, ST, 
396:      + CIK, CII, CKK, DT1, DT2, DT3, DT4, DT5, DT6, DT7, DT8, DT9,  STH409:      + CIK, CII, CKK, DT1, DT2, DT3, DT4, DT5, DT6, DT7, DT8, DT9,  STH
397: 410: 
398: 411: 
399:         DOUBLE PRECISION ANTC(NTA), Tk(NTA)412:         DOUBLE PRECISION ANTC(NTA), Tk(NTA)
400:         INTEGER II, IT(NTA), JT(NTA), KT(NTA)413:         INTEGER II, IT(NTA), JT(NTA), KT(NTA)
401:         INTEGER I3, J3, K3414:         INTEGER I3, J3, K3
402: 415: 
 416:       do II=1,NATOMS
 417:          gradOMP(II*3)=0
 418:          gradOMP(II*3-1)=0
 419:          gradOMP(II*3-2)=0
 420:       enddo
 421: 
 422:         etemp=0
403: !$OMP PARALLEL PRIVATE(RIJ,RKJ,RIK,ANT,XIJ,YIJ,ZIJ,XKJ,YKJ,423: !$OMP PARALLEL PRIVATE(RIJ,RKJ,RIK,ANT,XIJ,YIJ,ZIJ,XKJ,YKJ,
404: !$OMP&  ZKJ,DF,CT0,CT1,CT2,DA,ST,CIK,CII,CKK,DT1,DT2,DT3,DT4,424: !$OMP&  ZKJ,DF,CT0,CT1,CT2,DA,ST,CIK,CII,CKK,DT1,DT2,DT3,DT4,
405: !$OMP&  DT5,DT6,DT7,DT8,DT9,STH) REDUCTION(+:ENERGY,grad)425: !$OMP&  DT5,DT6,DT7,DT8,DT9,STH) REDUCTION(+:etemp,gradOMP)
406: !$OMP DO426: !$OMP DO
407:           DO II = 1, nTA427:           DO II = 1, nTA
408:             I3 = IT(II)428:             I3 = IT(II)
409:             J3 = JT(II)429:             J3 = JT(II)
410:             K3 = KT(II)430:             K3 = KT(II)
411:             XIJ = X(I3)-X(J3)431:             XIJ = X(I3)-X(J3)
412:             YIJ = Y(I3)-Y(J3)432:             YIJ = Y(I3)-Y(J3)
413:             ZIJ = Z(I3)-Z(J3)433:             ZIJ = Z(I3)-Z(J3)
414:             XKJ = X(K3)-X(J3)434:             XKJ = X(K3)-X(J3)
415:             YKJ = Y(K3)-Y(J3)435:             YKJ = Y(K3)-Y(J3)
430:             CKK = STH/RKJ450:             CKK = STH/RKJ
431:             DT1 = CIK*XKJ-CII*XIJ451:             DT1 = CIK*XKJ-CII*XIJ
432:             DT2 = CIK*YKJ-CII*YIJ452:             DT2 = CIK*YKJ-CII*YIJ
433:             DT3 = CIK*ZKJ-CII*ZIJ453:             DT3 = CIK*ZKJ-CII*ZIJ
434:             DT7 = CIK*XIJ-CKK*XKJ454:             DT7 = CIK*XIJ-CKK*XKJ
435:             DT8 = CIK*YIJ-CKK*YKJ455:             DT8 = CIK*YIJ-CKK*YKJ
436:             DT9 = CIK*ZIJ-CKK*ZKJ456:             DT9 = CIK*ZIJ-CKK*ZKJ
437:             DT4 = -DT1-DT7457:             DT4 = -DT1-DT7
438:             DT5 = -DT2-DT8458:             DT5 = -DT2-DT8
439:             DT6 = -DT3-DT9459:             DT6 = -DT3-DT9
440:             grad(I3*3-2) = grad(I3*3-2)+ DT1460:             gradOMP(I3*3-2) = gradOMP(I3*3-2)+ DT1
441:             grad(I3*3-1) = grad(I3*3-1)+ DT2461:             gradOMP(I3*3-1) = gradOMP(I3*3-1)+ DT2
442:             grad(I3*3)   = grad(I3*3)  + DT3462:             gradOMP(I3*3)   = gradOMP(I3*3)  + DT3
443:             grad(J3*3-2) = grad(J3*3-2)+ DT4463:             gradOMP(J3*3-2) = gradOMP(J3*3-2)+ DT4
444:             grad(J3*3-1) = grad(J3*3-1)+ DT5464:             gradOMP(J3*3-1) = gradOMP(J3*3-1)+ DT5
445:             grad(J3*3)   = grad(J3*3)  + DT6465:             gradOMP(J3*3)   = gradOMP(J3*3)  + DT6
446:             grad(K3*3-2) = grad(K3*3-2)+ DT7466:             gradOMP(K3*3-2) = gradOMP(K3*3-2)+ DT7
447:             grad(K3*3-1) = grad(K3*3-1)+ DT8467:             gradOMP(K3*3-1) = gradOMP(K3*3-1)+ DT8
448:             grad(K3*3)   = grad(K3*3)  + DT9468:             gradOMP(K3*3)   = gradOMP(K3*3)  + DT9
449:             ENERGY = ENERGY + TK(II)*(ANTC(II)- ANT)**2/2.0469:             etemp = etemp + TK(II)*(ANTC(II)- ANT)**2/2.0
450:           END DO470:           END DO
451: !$OMP END DO471: !$OMP END DO
452: !$OMP END PARALLEL472: !$OMP END PARALLEL
 473:         energy=energy+etemp 
 474:         do C1=1,NATOMS
 475:               grad(C1*3-2) = grad(C1*3-2)+ gradOMP(C1*3-2) 
 476:               grad(C1*3-1) = grad(C1*3-1)+ gradOMP(C1*3-1) 
 477:               grad(C1*3)   = grad(C1*3)  + gradOMP(C1*3)   
 478:         enddo
 479: 
453: 480: 
454:        RETURN481:        RETURN
455:        END482:        END
456: 483: 
457: !^^^^^^^^^^^^^^^^^^^^^^^^End of SBMANGL^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^484: !^^^^^^^^^^^^^^^^^^^^^^^^End of SBMANGL^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
458: 485: 
459: 486: 
460: 487: 
461: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<488: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
462: !* SBMdihedral computes the dihedral angles and the forces due to them *489: !* SBMdihedral computes the dihedral angles and the forces due to them *
463: !**********************************************************************490: !**********************************************************************
464: 491: 
465:       SUBROUTINE SBMdihedral(x,y,z,grad, energy, NATOMS,IP,JP,KP,LP,PK,492:       SUBROUTINE SBMdihedral(x,y,z,grad, energy, NATOMS,IP,JP,KP,LP,PK,
466:      Q PHITYPE,PHISBM,NPA)493:      Q PHITYPE,PHISBM,NPA)
467:       USE KEY494:       USE KEY
468:       implicit NONE495:       implicit NONE
469:       integer I, N, J, NATOMS, NPA, II496:       integer I, N, J, NATOMS, NPA, II
470:       DOUBLE PRECISION x(NATOMS),y(NATOMS),z(NATOMS),497:       DOUBLE PRECISION x(NATOMS),y(NATOMS),z(NATOMS),
471:      Q grad(3*NATOMS),energy498:      Q grad(3*NATOMS),gradOMP(3*NATOMS),energy,ETEMP
472: 499: 
473: 500: 
474:       DOUBLE PRECISION PK(NPA),PHISBM(NPA)501:       DOUBLE PRECISION PK(NPA),PHISBM(NPA)
475:       INTEGER IP(NPA), JP(NPA), KP(NPA), LP(NPA) 502:       INTEGER IP(NPA), JP(NPA), KP(NPA), LP(NPA) 
476:       INTEGER PHITYPE(NPA)503:       INTEGER PHITYPE(NPA)
477: 504: 
478:       double precision lfac505:       double precision lfac
479:       integer I3, J3, K3, L3,C1506:       integer I3, J3, K3, L3,C1
480:       double precision  XIJ,YIJ,ZIJ,XKJ,YKJ,507:       double precision  XIJ,YIJ,ZIJ,XKJ,YKJ,
481:      + ZKJ,XKL,YKL,ZKL,RIJ, RKJ,RKL,DX,DY,508:      + ZKJ,XKL,YKL,ZKL,RIJ, RKJ,RKL,DX,DY,
491:       double precision  TM24,TM06,tenm3,zero,one,two,four,six,twelve518:       double precision  TM24,TM06,tenm3,zero,one,two,four,six,twelve
492: 519: 
493:       DOUBLE PRECISION TT1, TT2, TT3, TT4, TT1X,TT1Y,TT1Z,TT2X,TT2Y,520:       DOUBLE PRECISION TT1, TT2, TT3, TT4, TT1X,TT1Y,TT1Z,TT2X,TT2Y,
494:      + TT2Z, TT3X, TT3Y, TT3Z, TT4X, TT4Y, TT4Z521:      + TT2Z, TT3X, TT3Y, TT3Z, TT4X, TT4Y, TT4Z
495: 522: 
496:       DATA TM24,TM06,tenm3/1.0d-24,1.0d-06,1.0d-03/523:       DATA TM24,TM06,tenm3/1.0d-24,1.0d-06,1.0d-03/
497:       data zero,one,two,four,six,twelve/0.d0,1.d0,2.d0,4.d0,6.d0,12.d0/524:       data zero,one,two,four,six,twelve/0.d0,1.d0,2.d0,4.d0,6.d0,12.d0/
498: 525: 
499:       double precision pi526:       double precision pi
500:       pi = 3.14159265358979323846264338327950288419716939937510527:       pi = 3.14159265358979323846264338327950288419716939937510
 528:       do i=1,NATOMS
 529:          gradOMP(i*3)=0
 530:          gradOMP(i*3-1)=0
 531:          gradOMP(i*3-2)=0
 532:       enddo
501: 533: 
 534:       etemp=0
502: !$OMP PARALLEL PRIVATE(I,I3,J3,K3,L3,II,XIJ,YIJ,ZIJ,XKJ,YKJ, 535: !$OMP PARALLEL PRIVATE(I,I3,J3,K3,L3,II,XIJ,YIJ,ZIJ,XKJ,YKJ, 
503: !$OMP& ZKJ,XKL,YKL,ZKL,RIJ,RKJ,RKL,DX,DY, 536: !$OMP& ZKJ,XKL,YKL,ZKL,RIJ,RKJ,RKL,DX,DY, 
504: !$OMP& DZ, GX,GY,GZ,CT,CPHI,537: !$OMP& DZ, GX,GY,GZ,CT,CPHI,
505: !$OMP& SPHI,Z1, Z2,FXI,FYI,FZI,538: !$OMP& SPHI,Z1, Z2,FXI,FYI,FZI,
506: !$OMP& FXJ,FYJ,FZJ, FXK,FYK,FZK,539: !$OMP& FXJ,FYJ,FZJ, FXK,FYK,FZK,
507: !$OMP& FXL,FYL,FZL,DF,Z10,Z20,Z12,Z11,Z22,CT0,CT1,AP0,AP1,540: !$OMP& FXL,FYL,FZL,DF,Z10,Z20,Z12,Z11,Z22,CT0,CT1,AP0,AP1,
508: !$OMP& Dums,DFLIM, DF1, DF0, DR1, DR2,DR3,DR4,DR5,DR6,DRX,DRY,DRZ,541: !$OMP& Dums,DFLIM, DF1, DF0, DR1, DR2,DR3,DR4,DR5,DR6,DRX,DRY,DRZ,
509: !$OMP& S,HGoverG,FGoverG,A1,A3,TT1,TT2,TT3,TT4,TT1X,TT1Y,TT1Z,TT2X,542: !$OMP& S,HGoverG,FGoverG,A1,A3,TT1,TT2,TT3,TT4,TT1X,TT1Y,TT1Z,TT2X,
510: !$OMP& TT2Y,TT2Z,TT3X,TT3Y,TT3Z,TT4X,TT4Y,TT4Z) REDUCTION(+:energy,grad)543: !$OMP& TT2Y,TT2Z,TT3X,TT3Y,TT3Z,TT4X,TT4Y,TT4Z) REDUCTION(+:etemp,gradOMP)
511: !!544: !!
512: !$OMP DO545: !$OMP DO
513:           DO II = 1,nPA546:           DO II = 1,nPA
514: 547: 
515:             I3 = IP(II)548:             I3 = IP(II)
516:             J3 = JP(II)549:             J3 = JP(II)
517:             K3 = KP(II)550:             K3 = KP(II)
518:             L3 = LP(II)551:             L3 = LP(II)
519:             XIJ = X(I3)-X(J3)552:             XIJ = X(I3)-X(J3)
520:             YIJ = Y(I3)-Y(J3)553:             YIJ = Y(I3)-Y(J3)
549:             AP0 = ACOS(CT1)582:             AP0 = ACOS(CT1)
550:             AP1 = PI-SIGN(AP0,S)583:             AP1 = PI-SIGN(AP0,S)
551:             CT = AP1584:             CT = AP1
552:             CPHI = COS(AP1)585:             CPHI = COS(AP1)
553:             SPHI = SIN(AP1)586:             SPHI = SIN(AP1)
554: ! Here is the energy part587: ! Here is the energy part
555:           A1=CT-PHISBM(II)588:           A1=CT-PHISBM(II)
556:           A3=A1*3589:           A3=A1*3
557: 590: 
558:         if(PHITYPE(II) .eq. 1)then591:         if(PHITYPE(II) .eq. 1)then
559:           energy =  energy + PK(II)*(3.0/2.0-cos(A1)-0.5*cos(A3))592:           etemp =  etemp + PK(II)*(3.0/2.0-cos(A1)-0.5*cos(A3))
560: ! dE/dPHI593: ! dE/dPHI
561:             DF=PK(II)*(sin(A1)+1.5*sin(A3))594:             DF=PK(II)*(sin(A1)+1.5*sin(A3))
562:         elseif(PHITYPE(II) .eq. 2)then595:         elseif(PHITYPE(II) .eq. 2)then
563:           if(A1 .gt. PI)then596:           if(A1 .gt. PI)then
564:                 A1=A1-2*PI597:                 A1=A1-2*PI
565:           elseif(A1 .lt. -PI)then598:           elseif(A1 .lt. -PI)then
566:                 A1=A1+2*PI599:                 A1=A1+2*PI
567:           endif600:           endif
568: 601: 
569:           energy =  energy + 0.5*PK(II)*A1**2602:           etemp =  etemp + 0.5*PK(II)*A1**2
570: ! dE/dPHI603: ! dE/dPHI
571:             DF=PK(II)*A1604:             DF=PK(II)*A1
572: 605: 
573:         else606:         else
574: 607: 
575:         write(*,*) 'unrecognized dihedral type', PHITYPE(II)608:         write(*,*) 'unrecognized dihedral type', PHITYPE(II)
576:         STOP609:         STOP
577:         endif610:         endif
578: 611: 
579: ! insert the new 612: ! insert the new 
594:           TT1Z=TT1*DZ627:           TT1Z=TT1*DZ
595:           TT2X=TT2*DX628:           TT2X=TT2*DX
596:           TT2Y=TT2*DY629:           TT2Y=TT2*DY
597:           TT2Z=TT2*DZ630:           TT2Z=TT2*DZ
598:           TT3X=TT3*GX631:           TT3X=TT3*GX
599:           TT3Y=TT3*GY632:           TT3Y=TT3*GY
600:           TT3Z=TT3*GZ633:           TT3Z=TT3*GZ
601:           TT4X=TT4*GX634:           TT4X=TT4*GX
602:           TT4Y=TT4*GY635:           TT4Y=TT4*GY
603:           TT4Z=TT4*GZ636:           TT4Z=TT4*GZ
604:           grad(I3*3-2) =  grad(I3*3-2)  + TT1X  637:           gradOMP(I3*3-2) =  gradOMP(I3*3-2)  + TT1X  
605:           grad(I3*3-1) =  grad(I3*3-1)  + TT1Y638:           gradOMP(I3*3-1) =  gradOMP(I3*3-1)  + TT1Y
606:           grad(I3*3)   =  grad(I3*3)    + TT1Z639:           gradOMP(I3*3)   =  gradOMP(I3*3)    + TT1Z
607:           grad(J3*3-2) =  grad(J3*3-2)  - TT1X - TT2X - TT3X640:           gradOMP(J3*3-2) =  gradOMP(J3*3-2)  - TT1X - TT2X - TT3X
608:           grad(J3*3-1) =  grad(J3*3-1)  - TT1Y - TT2Y - TT3Y641:           gradOMP(J3*3-1) =  gradOMP(J3*3-1)  - TT1Y - TT2Y - TT3Y
609:           grad(J3*3)   =  grad(J3*3)    - TT1Z - TT2Z - TT3Z642:           gradOMP(J3*3)   =  gradOMP(J3*3)    - TT1Z - TT2Z - TT3Z
610:           grad(K3*3-2) =  grad(K3*3-2)  + TT2X + TT3X - TT4X643:           gradOMP(K3*3-2) =  gradOMP(K3*3-2)  + TT2X + TT3X - TT4X
611:           grad(K3*3-1) =  grad(K3*3-1)  + TT2Y + TT3Y - TT4Y644:           gradOMP(K3*3-1) =  gradOMP(K3*3-1)  + TT2Y + TT3Y - TT4Y
612:           grad(K3*3)   =  grad(K3*3)    + TT2Z + TT3Z - TT4Z645:           gradOMP(K3*3)   =  gradOMP(K3*3)    + TT2Z + TT3Z - TT4Z
613:           grad(L3*3-2) =  grad(L3*3-2)  + TT4X646:           gradOMP(L3*3-2) =  gradOMP(L3*3-2)  + TT4X
614:           grad(L3*3-1) =  grad(L3*3-1)  + TT4Y647:           gradOMP(L3*3-1) =  gradOMP(L3*3-1)  + TT4Y
615:           grad(L3*3)   =  grad(L3*3)    + TT4Z648:           gradOMP(L3*3)   =  gradOMP(L3*3)    + TT4Z
616:         END DO649:         END DO
617: !$OMP END DO 650: !$OMP END DO 
618: !$OMP END PARALLEL651: !$OMP END PARALLEL
619: 652: 
 653:         ENERGY=ENERGY+etemp
 654:         do C1=1,NATOMS
 655:               grad(C1*3-2) = grad(C1*3-2)+ gradOMP(C1*3-2) 
 656:               grad(C1*3-1) = grad(C1*3-1)+ gradOMP(C1*3-1) 
 657:               grad(C1*3)   = grad(C1*3)  + gradOMP(C1*3)   
 658:         enddo
 659: 
 660: 
620:           END661:           END
621: 662: 
622: 663: 
623: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^END of SBMDihedral^^^^^^^^^^^^^^^^^^^^^^^^^^^664: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^END of SBMDihedral^^^^^^^^^^^^^^^^^^^^^^^^^^^
624: 665: 
625: 666: 
626: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<667: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
627: !* SBMCONTACTS: computes the force on all atoms due to contacts via a   *668: !* SBMCONTACTS: computes the force on all atoms due to contacts via a   *
628: !* 10-12 or 6-12 potential                                              *669: !* 10-12 or 6-12 potential                                              *
629: !***********************************************************************670: !***********************************************************************
630: 671: 
631:       subroutine SBMcontacts(x,y,z,grad,energy,672:       subroutine SBMcontacts(x,y,z,grad,energy,
632:      Q NATOMS,IC,JC,Sigma,EpsC,NC,CONTACTTYPE)673:      Q NATOMS,IC,JC,Sigma,EpsC,NC,CONTACTTYPE)
633:       USE KEY674:       USE KEY
634:       implicit NONE675:       implicit NONE
635:       integer I, NATOMS,NC676:       integer I, NATOMS,NC
636: 677: 
637:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS), 678:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS), 
638:      Q grad(3*NATOMS), energy679:      Q grad(3*NATOMS), energy,etemp, GRADOMP(3*NATOMS)
639:       DOUBLE PRECISION dx,dy,dz680:       DOUBLE PRECISION dx,dy,dz
640: 681: 
641:       integer C1, C2 682:       integer C1, C2 
642:       DOUBLE PRECISION  r2, rm2, rm10, rm6,f_over_r,r1 683:       DOUBLE PRECISION  r2, rm2, rm10, rm6,f_over_r,r1 
643: 684: 
644:         DOUBLE PRECISION Sigma(NC), EpsC(NC)685:         DOUBLE PRECISION Sigma(NC), EpsC(NC)
645:         INTEGER IC(NC), JC(NC),CONTACTTYPE686:         INTEGER IC(NC), JC(NC),CONTACTTYPE
646: 687: 
647: ! type 1 is 6-12688: ! type 1 is 6-12
648: ! type 2 is 12-12689: ! type 2 is 12-12
649: ! implement type 3 as gaussian690: ! implement type 3 as gaussian
650: ! implement type 4 as dual-basin gaussian691: ! implement type 4 as dual-basin gaussian
651: ! question: Should type be part of the interaction, or should they be organized692: ! question: Should type be part of the interaction, or should they be organized
652: ! into separate lists?693: ! into separate lists?
 694:       do I=1,NATOMS
 695:          gradOMP(I*3)=0
 696:          gradOMP(I*3-1)=0
 697:          gradOMP(I*3-2)=0
 698:       enddo
 699: 
653: 700: 
 701:         ETEMP=0
654: ! type 1 is 6-12 interaction702: ! type 1 is 6-12 interaction
655:         if(CONTACTTYPE .eq. 1)then703:         if(CONTACTTYPE .eq. 1)then
656: !$OMP PARALLEL PRIVATE(I,C1,C2,DX,DY,DZ,R2,RM2,RM6,F_OVER_R)REDUCTION(+:ENERGY,grad)704: !$OMP PARALLEL PRIVATE(I,C1,C2,DX,DY,DZ,R2,RM2,RM6,F_OVER_R)REDUCTION(+:etemp,gradOMP)
657: !$OMP DO705: !$OMP DO
658:           do i=1, NC706:           do i=1, NC
659:             C1 = IC(i)707:             C1 = IC(i)
660:             C2 = JC(i)708:             C2 = JC(i)
661:             dx = X(C1) - X(C2)709:             dx = X(C1) - X(C2)
662:             dy = Y(C1) - Y(C2)710:             dy = Y(C1) - Y(C2)
663:             dz = Z(C1) - Z(C2)711:             dz = Z(C1) - Z(C2)
664:             r2 = dx**2 + dy**2 + dz**2712:             r2 = dx**2 + dy**2 + dz**2
665:             rm2 = 1.0/r2713:             rm2 = 1.0/r2
666:             rm2 = rm2*sigma(i)714:             rm2 = rm2*sigma(i)
667:             rm6 = rm2**3715:             rm6 = rm2**3
668:             ENERGY = ENERGY + epsC(i)*rm6*(rm6-2.0)716:             ETEMP = ETEMP + epsC(i)*rm6*(rm6-2.0)
669:             f_over_r = -epsC(i)*12.0*rm6*(rm6-1.0)/r2717:             f_over_r = -epsC(i)*12.0*rm6*(rm6-1.0)/r2
670:             grad(3*C1-2) = grad(3*C1-2) + f_over_r * dx718:             gradOMP(3*C1-2) = gradOMP(3*C1-2) + f_over_r * dx
671:             grad(3*C1-1) = grad(3*C1-1) + f_over_r * dy719:             gradOMP(3*C1-1) = gradOMP(3*C1-1) + f_over_r * dy
672:             grad(3*C1)   = grad(3*C1)   + f_over_r * dz720:             gradOMP(3*C1)   = gradOMP(3*C1)   + f_over_r * dz
673:             grad(3*C2-2) = grad(3*C2-2) - f_over_r * dx721:             gradOMP(3*C2-2) = gradOMP(3*C2-2) - f_over_r * dx
674:             grad(3*C2-1) = grad(3*C2-1) - f_over_r * dy722:             gradOMP(3*C2-1) = gradOMP(3*C2-1) - f_over_r * dy
675:             grad(3*C2)   = grad(3*C2)   - f_over_r * dz723:             gradOMP(3*C2)   = gradOMP(3*C2)   - f_over_r * dz
676:           enddo724:           enddo
677: !$OMP END DO725: !$OMP END DO
678: !$OMP END PARALLEL726: !$OMP END PARALLEL
679: ! type 2 is 10-12 interaction727: ! type 2 is 10-12 interaction
680:         elseif(CONTACTTYPE .eq. 2)then728:         elseif(CONTACTTYPE .eq. 2)then
681: !$OMP PARALLEL PRIVATE(I,C1,C2,DX,DY,DZ,R2,RM2,RM10,F_OVER_R)REDUCTION(+:ENERGY,grad)729: !$OMP PARALLEL PRIVATE(I,C1,C2,DX,DY,DZ,R2,RM2,RM10,F_OVER_R)REDUCTION(+:etemp,gradOMP)
682: !$OMP DO730: !$OMP DO
683:           do i=1, NC731:           do i=1, NC
684:              C1 = IC(i)732:              C1 = IC(i)
685:              C2 = JC(i)733:              C2 = JC(i)
686:              dx = X(C1) - X(C2)734:              dx = X(C1) - X(C2)
687:              dy = Y(C1) - Y(C2)735:              dy = Y(C1) - Y(C2)
688:              dz = Z(C1) - Z(C2)736:              dz = Z(C1) - Z(C2)
689:              r2 = dx**2 + dy**2 + dz**2737:              r2 = dx**2 + dy**2 + dz**2
690:              rm2 = 1.0/r2738:              rm2 = 1.0/r2
691:              rm2 = rm2*sigma(i)739:              rm2 = rm2*sigma(i)
692:              rm10 = rm2**5740:              rm10 = rm2**5
693:              ENERGY = ENERGY + epsC(i)*rm10*(5*rm2-6.0)741:              ETEMP = ETEMP + epsC(i)*rm10*(5*rm2-6.0)
694:              f_over_r = -epsc(i)*60.0*rm10*(rm2-1.0)/r2742:              f_over_r = -epsc(i)*60.0*rm10*(rm2-1.0)/r2
695:              grad(3*C1-2) = grad(3*C1-2) + f_over_r * dx743:              gradOMP(3*C1-2) = gradOMP(3*C1-2) + f_over_r * dx
696:              grad(3*C1-1) = grad(3*C1-1) + f_over_r * dy744:              gradOMP(3*C1-1) = gradOMP(3*C1-1) + f_over_r * dy
697:              grad(3*C1)   = grad(3*C1)   + f_over_r * dz745:              gradOMP(3*C1)   = gradOMP(3*C1)   + f_over_r * dz
698:              grad(3*C2-2) = grad(3*C2-2) - f_over_r * dx746:              gradOMP(3*C2-2) = gradOMP(3*C2-2) - f_over_r * dx
699:              grad(3*C2-1) = grad(3*C2-1) - f_over_r * dy747:              gradOMP(3*C2-1) = gradOMP(3*C2-1) - f_over_r * dy
700:              grad(3*C2)   = grad(3*C2)   - f_over_r * dz748:              gradOMP(3*C2)   = gradOMP(3*C2)   - f_over_r * dz
701:            enddo749:            enddo
702: !$OMP END DO750: !$OMP END DO
703: !$OMP END PARALLEL751: !$OMP END PARALLEL
704: 752: 
705:         else753:         else
706:         write(*,*) CONTACTTYPE, 'is not a valid contact selection'754:         write(*,*) CONTACTTYPE, 'is not a valid contact selection'
707:         stop755:         stop
708:         endif756:         endif
709: 757: 
 758:         ENERGY=ENERGY+ETEMP
 759:         do I=1,NATOMS
 760:               grad(I*3-2) = grad(I*3-2)+ gradOMP(I*3-2) 
 761:               grad(I*3-1) = grad(I*3-1)+ gradOMP(I*3-1) 
 762:               grad(I*3)   = grad(I*3)  + gradOMP(I*3)   
 763:         enddo
 764: 
 765: 
710: 766: 
711:       end767:       end
712: 768: 
713: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^end of SBMContacts^^^^^^^^^^^^^^^^^^^^^^^^^^^769: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^end of SBMContacts^^^^^^^^^^^^^^^^^^^^^^^^^^^
714: 770: 
715: 771: 
716: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<772: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
717: !* SBMNonContacts computes the forces due to non native contacts       *773: !* SBMNonContacts computes the forces due to non native contacts       *
718: !**********************************************************************774: !**********************************************************************
719: 775: 
720:       subroutine SBMnoncontacts(x,y,z,grad, energy, 776:       subroutine SBMnoncontacts(x,y,z,grad, energy, 
721:      Q NATOMS,NNCinc,NCswitch,NCcut,STT)777:      Q NATOMS,NNCinc,NCswitch,NCcut,STT)
722:       USE KEY778:       USE KEY
723:       implicit NONE779:       implicit NONE
724:       integer I, N, J, AN, NATOMS780:       integer I, N, J, AN, NATOMS
725: 781: 
726: 782: 
727:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS),783:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS),
728:      Q grad(3*NATOMS), energy,STT,SA,SB,SC784:      Q grad(3*NATOMS), gradOMP(3*NATOMS), energy,STT,SA,SB,SC
729: 785: 
730:       integer C1, C2, ii,jj,kk,k,l,iii,jjj786:       integer C1, C2, ii,jj,kk,k,l,iii,jjj
731:       DOUBLE PRECISION  r2, rm2, rm14, f_over_r, NCswitch,NCcut 787:       DOUBLE PRECISION  r2, rm2, rm14, f_over_r, NCswitch,NCcut 
732:         INTEGER NTHREADS788:         INTEGER NTHREADS
733: !$    INTEGER  OMP_GET_NUM_THREADS,789: !$    INTEGER  OMP_GET_NUM_THREADS,
734: !$   Q  OMP_GET_THREAD_NUM,OMP_NUM_THREADS790: !$   Q  OMP_GET_THREAD_NUM,OMP_NUM_THREADS
735: 791: 
736:         INTEGER NNCinc(NATOMS,NATOMS)792:         INTEGER NNCinc(NATOMS,NATOMS)
737:         DOUBLE PRECISION dx,dy,dz793:         DOUBLE PRECISION dx,dy,dz
738:         integer tempN, alpha794:         integer tempN, alpha
744:         ! number of atoms per grid, max800:         ! number of atoms per grid, max
745:         parameter (maxpergrid=50)801:         parameter (maxpergrid=50)
746:         ! dimensions of grid802:         ! dimensions of grid
747:         parameter (maxgrid=15)803:         parameter (maxgrid=15)
748:         dimension Ngrid(maxgrid,maxgrid,maxgrid),804:         dimension Ngrid(maxgrid,maxgrid,maxgrid),
749:      Q  grid(maxgrid,maxgrid,maxgrid,maxpergrid)805:      Q  grid(maxgrid,maxgrid,maxgrid,maxpergrid)
750:         integer MaxGridX,MaxGridY,MaxGridZ806:         integer MaxGridX,MaxGridY,MaxGridZ
751:         double precision gridsize,RD1807:         double precision gridsize,RD1
752:         double precision minX,minY,minZ,maxX,maxY,maxZ808:         double precision minX,minY,minZ,maxX,maxY,maxZ
753:         integer Xgrid,Ygrid,Zgrid,TID809:         integer Xgrid,Ygrid,Zgrid,TID
 810:         double precision Etemp
754:         Rdiff=NCcut-NCswitch811:         Rdiff=NCcut-NCswitch
755:         alpha=12812:         alpha=12
756: 813: 
757: 814: 
758:         GRIDSIZE=NCcut*1.01815:         GRIDSIZE=NCcut*1.01
759:         Rcut2=NCcut**2816:         Rcut2=NCcut**2
760:         Rswitch2=NCswitch**2817:         Rswitch2=NCswitch**2
761: 818: 
762:         SB=-1.0/Rdiff**3*( 2*alpha*STT/NCcut**(alpha+1)  + 819:         SB=-1.0/Rdiff**3*( 2*alpha*STT/NCcut**(alpha+1)  + 
763:      Q    (alpha)*(alpha+1)*STT*Rdiff/NCcut**(alpha+2)   )820:      Q    (alpha)*(alpha+1)*STT*Rdiff/NCcut**(alpha+2)   )
770:         minX=10000000827:         minX=10000000
771:         minY=10000000828:         minY=10000000
772:         minZ=10000000829:         minZ=10000000
773:         830:         
774:         maxX=-10000000831:         maxX=-10000000
775:         maxY=-10000000832:         maxY=-10000000
776:         maxZ=-10000000833:         maxZ=-10000000
777: 834: 
778:         do i=1,NATOMS835:         do i=1,NATOMS
779: 836: 
 837:            gradOMP(i*3)=0
 838:            gradOMP(i*3-1)=0
 839:            gradOMP(i*3-2)=0
780:            if(X(i) .gt. maxX)then840:            if(X(i) .gt. maxX)then
781:                 maxX=X(i)841:                 maxX=X(i)
782:            elseif(Y(i) .gt. maxY)then842:            elseif(Y(i) .gt. maxY)then
783:                 maxY=Y(i)843:                 maxY=Y(i)
784:            endif844:            endif
785: 845: 
786:            if(Z(i) .gt. maxZ)then846:            if(Z(i) .gt. maxZ)then
787:                 maxZ=Z(i)847:                 maxZ=Z(i)
788:            elseif(X(i) .lt. minX)then848:            elseif(X(i) .lt. minX)then
789:                 minX=X(i)849:                 minX=X(i)
831:                 grid(Xgrid,Ygrid,Zgrid,Ngrid(Xgrid,Ygrid,Zgrid))=i891:                 grid(Xgrid,Ygrid,Zgrid,Ngrid(Xgrid,Ygrid,Zgrid))=i
832:         enddo892:         enddo
833: 893: 
834: 894: 
835:            tempN=0895:            tempN=0
836: 896: 
837: ! add a second loop that goes over charged atoms897: ! add a second loop that goes over charged atoms
838: 898: 
839: !$OMP PARALLEL899: !$OMP PARALLEL
840: !$OMP& PRIVATE(i,ii,jj,kk,jjj,C1,C2,r2,rm2,rm14,f_over_r,RD1,dx,dy,dz,Xgrid,Ygrid,Zgrid,TID)  900: !$OMP& PRIVATE(i,ii,jj,kk,jjj,C1,C2,r2,rm2,rm14,f_over_r,RD1,dx,dy,dz,Xgrid,Ygrid,Zgrid,TID)  
841: !$OMP& REDUCTION(+:ENERGY,grad)901: !$OMP& REDUCTION(+:etemp,gradOMP)
 902:         etemp=0
842:         TID=0903:         TID=0
843:         NTHREADS=1;904:         NTHREADS=1;
844: !$      TID = OMP_GET_THREAD_NUM()905: !$      TID = OMP_GET_THREAD_NUM()
845: !$      NTHREADS = OMP_GET_NUM_THREADS()906: !$      NTHREADS = OMP_GET_NUM_THREADS()
846:         do i=1+TID,NATOMS,NTHREADS907:         do i=1+TID,NATOMS,NTHREADS
847: 908: 
848: !!!        do i=1,NATOMS909: !!!        do i=1,NATOMS
849:            C1 = i910:            C1 = i
850:                 Xgrid=int((X(i)-minX)/gridsize)+1911:                 Xgrid=int((X(i)-minX)/gridsize)+1
851:                 Ygrid=int((Y(i)-minY)/gridsize)+1912:                 Ygrid=int((Y(i)-minY)/gridsize)+1
862:            C2=grid(ii,jj,kk,jjj)923:            C2=grid(ii,jj,kk,jjj)
863:              if(NNCinc(i,C2) .eq. 1)then924:              if(NNCinc(i,C2) .eq. 1)then
864: !!! SSS925: !!! SSS
865:         dx = X(C1) - X(C2)926:         dx = X(C1) - X(C2)
866:         dy = Y(C1) - Y(C2)927:         dy = Y(C1) - Y(C2)
867:         dz = Z(C1) - Z(C2)928:         dz = Z(C1) - Z(C2)
868:           r2 = dx**2 + dy**2 + dz**2929:           r2 = dx**2 + dy**2 + dz**2
869:           if(r2 .le. Rcut2)then930:           if(r2 .le. Rcut2)then
870:              rm2 = 1/r2931:              rm2 = 1/r2
871:              rm14 = rm2**7932:              rm14 = rm2**7
872:                 energy=energy+STT*rm2**6+SC933:                 etemp=etemp+STT*rm2**6+SC
873:                 f_over_r=-STT*12.0*rm14934:                 f_over_r=-STT*12.0*rm14
874:                 RD1=sqrt(r2)-NCswitch935:                 RD1=sqrt(r2)-NCswitch
875:                 if(r2 .gt. Rswitch2)then936:                 if(r2 .gt. Rswitch2)then
876:                         f_over_r=f_over_r+(SA*RD1**2+SB*RD1**3)*sqrt(rm2)937:                         f_over_r=f_over_r+(SA*RD1**2+SB*RD1**3)*sqrt(rm2)
877:                         energy=energy+SA/3.0*RD1**3+SB/4.0*RD1**4938:                         etemp=etemp+SA/3.0*RD1**3+SB/4.0*RD1**4
878:                 elseif(r2 .lt. Rswitch2)then939:                 elseif(r2 .lt. Rswitch2)then
879:                 ! normal repulsive term940:                 ! normal repulsive term
880:                 else941:                 else
881:                 ! things should have fallen in one of the previous two...942:                 ! things should have fallen in one of the previous two...
882:                 write(*,*) 'something went wrong with switching function'943:                 write(*,*) 'something went wrong with switching function'
883: !                call abort944: !                call abort
884:                 endif945:                 endif
885:               grad(C1*3-2) = grad(C1*3-2) + f_over_r * dx946:               gradOMP(C1*3-2) = gradOMP(C1*3-2) + f_over_r * dx
886:               grad(C1*3-1) = grad(C1*3-1) + f_over_r * dy947:               gradOMP(C1*3-1) = gradOMP(C1*3-1) + f_over_r * dy
887:               grad(C1*3)   = grad(C1*3)   + f_over_r * dz948:               gradOMP(C1*3)   = gradOMP(C1*3)   + f_over_r * dz
888:               grad(C2*3-2) = grad(C2*3-2) - f_over_r * dx949: 
889:               grad(C2*3-1) = grad(C2*3-1) - f_over_r * dy950:               gradOMP(C2*3-2) =  gradOMP(C2*3-2) - f_over_r * dx
890:               grad(C2*3)   = grad(C2*3)   - f_over_r * dz951:               gradOMP(C2*3-1) =  gradOMP(C2*3-1) - f_over_r * dy
 952:               gradOMP(C2*3)   =  gradOMP(C2*3)   - f_over_r * dz
891: 953: 
892:             endif954:             endif
893: 955: 
894: 956: 
895:            endif957:            endif
896:            enddo958:            enddo
897:  959:  
898:           endif960:           endif
899: 961: 
900:             enddo962:             enddo
901:            enddo963:            enddo
902:           enddo964:           enddo
903: 965: 
904:           enddo966:           enddo
905: !$OMP END PARALLEL967: !$OMP END PARALLEL
 968:           energy = energy  + etemp
 969:         do C1=1,NATOMS
 970:               grad(C1*3-2) = grad(C1*3-2)+ gradOMP(C1*3-2) 
 971:               grad(C1*3-1) = grad(C1*3-1)+ gradOMP(C1*3-1) 
 972:               grad(C1*3)   = grad(C1*3)  + gradOMP(C1*3)   
 973:         enddo
906: 974: 
907:       END975:       END
908: 976: 
909: !^^^^^^^^^^^^^^^^^^^^^^^^^^^End of SBMNonContacts^^^^^^^^^^^^^^^^^^^^^^^^^977: !^^^^^^^^^^^^^^^^^^^^^^^^^^^End of SBMNonContacts^^^^^^^^^^^^^^^^^^^^^^^^^
910: 978: 
911: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<979: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
912: ! FUNCTIONS for DH Calculations980: ! FUNCTIONS for DH Calculations
913: !*****************************************************981: !*****************************************************
914: 982: 
915:       DOUBLE PRECISION FUNCTION SBMDHV(kappa,r)983:       DOUBLE PRECISION FUNCTION SBMDHV(kappa,r)
933: 1001: 
934:       subroutine SBMDHELEC(x,y,z,grad, energy, natoms, 1002:       subroutine SBMDHELEC(x,y,z,grad, energy, natoms, 
935:      Q NNCinc,PREFACTOR,KAPPA,DHswitch,DHcut,SBMCHARGE,SBMCHARGEON)1003:      Q NNCinc,PREFACTOR,KAPPA,DHswitch,DHcut,SBMCHARGE,SBMCHARGEON)
936: 1004: 
937:       USE KEY1005:       USE KEY
938:       implicit NONE1006:       implicit NONE
939:       integer I, N, J, AN, NATOMS1007:       integer I, N, J, AN, NATOMS
940: 1008: 
941: 1009: 
942:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS),1010:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS),
943:      Q grad(3*NATOMS), energy,STT,SA,SB,SC,1011:      Q grad(3*NATOMS), gradOMP(3*NATOMS), energy,STT,SA,SB,SC,
944:      Q SBMDHVP, SBMDHV1012:      Q SBMDHVP, SBMDHV
945:       integer C1, C2, ii,jj,kk,k,l,iii,jjj1013:       integer C1, C2, ii,jj,kk,k,l,iii,jjj
946:       DOUBLE PRECISION  r1,r2, rm2, rm14, f_over_r, NCswitch,NCcut 1014:       DOUBLE PRECISION  r1,r2, rm2, rm14, f_over_r, NCswitch,NCcut 
947:       DOUBLE PRECISION A,B,C, D, COEF2, COEF31015:       DOUBLE PRECISION A,B,C, D, COEF2, COEF3
948:         INTEGER NTHREADS,  OMP_GET_NUM_THREADS,1016:         INTEGER NTHREADS,  OMP_GET_NUM_THREADS,
949:      Q  OMP_GET_THREAD_NUM,OMP_NUM_THREADS1017:      Q  OMP_GET_THREAD_NUM,OMP_NUM_THREADS
950: 1018: 
951:         INTEGER NNCinc(NATOMS,NATOMS)1019:         INTEGER NNCinc(NATOMS,NATOMS)
952:         DOUBLE PRECISION dx,dy,dz1020:         DOUBLE PRECISION dx,dy,dz
953:         integer tempN, alpha,SBMCHARGEON(NATOMS+1)1021:         integer tempN, alpha,SBMCHARGEON(NATOMS+1)
959:         ! number of atoms per grid, max1027:         ! number of atoms per grid, max
960:         parameter (maxpergrid=50)1028:         parameter (maxpergrid=50)
961:         ! dimensions of grid1029:         ! dimensions of grid
962:         parameter (maxgrid=15)1030:         parameter (maxgrid=15)
963:         dimension Ngrid(maxgrid,maxgrid,maxgrid),1031:         dimension Ngrid(maxgrid,maxgrid,maxgrid),
964:      Q  grid(maxgrid,maxgrid,maxgrid,maxpergrid)1032:      Q  grid(maxgrid,maxgrid,maxgrid,maxpergrid)
965:         integer MaxGridX,MaxGridY,MaxGridZ1033:         integer MaxGridX,MaxGridY,MaxGridZ
966:         double precision gridsize,RD11034:         double precision gridsize,RD1
967:         double precision minX,minY,minZ,maxX,maxY,maxZ1035:         double precision minX,minY,minZ,maxX,maxY,maxZ
968:         integer Xgrid,Ygrid,Zgrid,TID1036:         integer Xgrid,Ygrid,Zgrid,TID
969:         double precision C2T1037:         double precision Etemp,C2T
970: 1038: 
971: 1039: 
972:       if(SBMCHARGEON(1) .eq. 1) RETURN1040:       if(SBMCHARGEON(1) .eq. 1) RETURN
973: 1041: 
974:         diff=DHcut-DHswitch1042:         diff=DHcut-DHswitch
975:         alpha=121043:         alpha=12
976: 1044: 
977: 1045: 
978:         GRIDSIZE=DHCUT*1.011046:         GRIDSIZE=DHCUT*1.01
979:         DHcut2=DHcut**21047:         DHcut2=DHcut**2
990:         minX=100000001058:         minX=10000000
991:         minY=100000001059:         minY=10000000
992:         minZ=100000001060:         minZ=10000000
993:         1061:         
994:         maxX=-100000001062:         maxX=-10000000
995:         maxY=-100000001063:         maxY=-10000000
996:         maxZ=-100000001064:         maxZ=-10000000
997: 1065: 
998:         do ii=2,SBMCHARGEON(1)1066:         do ii=2,SBMCHARGEON(1)
999:            i=SBMCHARGEON(ii)1067:            i=SBMCHARGEON(ii)
 1068:            gradOMP(i*3)=0
 1069:            gradOMP(i*3-1)=0
 1070:            gradOMP(i*3-2)=0
1000: 1071: 
1001:            if(X(i) .gt. maxX)then1072:            if(X(i) .gt. maxX)then
1002:                 maxX=X(i)1073:                 maxX=X(i)
1003:            elseif(Y(i) .gt. maxY)then1074:            elseif(Y(i) .gt. maxY)then
1004:                 maxY=Y(i)1075:                 maxY=Y(i)
1005:            endif1076:            endif
1006: 1077: 
1007:            if(Z(i) .gt. maxZ)then1078:            if(Z(i) .gt. maxZ)then
1008:                 maxZ=Z(i)1079:                 maxZ=Z(i)
1009:            elseif(X(i) .lt. minX)then1080:            elseif(X(i) .lt. minX)then
1054:                 endif1125:                 endif
1055:                 grid(Xgrid,Ygrid,Zgrid,Ngrid(Xgrid,Ygrid,Zgrid))=i1126:                 grid(Xgrid,Ygrid,Zgrid,Ngrid(Xgrid,Ygrid,Zgrid))=i
1056:         enddo1127:         enddo
1057: 1128: 
1058: 1129: 
1059:            tempN=01130:            tempN=0
1060: 1131: 
1061: 1132: 
1062: !$OMP PARALLEL1133: !$OMP PARALLEL
1063: !$OMP& PRIVATE(i,ii,iii,jj,kk,jjj,C1,C2,r2,rm2,rm14,f_over_r,RD1,dx,dy,dz,Xgrid,Ygrid,Zgrid,TID,C2T)  1134: !$OMP& PRIVATE(i,ii,iii,jj,kk,jjj,C1,C2,r2,rm2,rm14,f_over_r,RD1,dx,dy,dz,Xgrid,Ygrid,Zgrid,TID,C2T)  
1064: !$OMP& REDUCTION(+:energy,grad)1135: !$OMP& REDUCTION(+:etemp,gradOMP)
 1136:         etemp=0
 1137:         etemp=0
1065:         TID=01138:         TID=0
1066:         NTHREADS=1;1139:         NTHREADS=1;
1067: !$      TID = OMP_GET_THREAD_NUM()1140: !$      TID = OMP_GET_THREAD_NUM()
1068: !$      NTHREADS = OMP_GET_NUM_THREADS()1141: !$      NTHREADS = OMP_GET_NUM_THREADS()
1069:         do iii=2+TID,SBMCHARGEON(1),NTHREADS1142:         do iii=2+TID,SBMCHARGEON(1),NTHREADS
1070: 1143: 
1071: 1144: 
1072:            C1=SBMCHARGEON(iii)1145:            C1=SBMCHARGEON(iii)
1073:                 Xgrid=int((X(C1)-minX)/gridsize)+11146:                 Xgrid=int((X(C1)-minX)/gridsize)+1
1074:                 Ygrid=int((Y(C1)-minY)/gridsize)+11147:                 Ygrid=int((Y(C1)-minY)/gridsize)+1
1088:         dx = X(C1) - X(C2)1161:         dx = X(C1) - X(C2)
1089:         dy = Y(C1) - Y(C2)1162:         dy = Y(C1) - Y(C2)
1090:         dz = Z(C1) - Z(C2)1163:         dz = Z(C1) - Z(C2)
1091:           r2 = dx**2 + dy**2 + dz**21164:           r2 = dx**2 + dy**2 + dz**2
1092:           if(r2 .le. DHcut2)then1165:           if(r2 .le. DHcut2)then
1093:            C2T=SBMCHARGE(C1)*SBMCHARGE(C2)1166:            C2T=SBMCHARGE(C1)*SBMCHARGE(C2)
1094:         ! add force evaluation now1167:         ! add force evaluation now
1095: 1168: 
1096:              rm2 = 1/r21169:              rm2 = 1/r2
1097:              r1=sqrt(r2)1170:              r1=sqrt(r2)
1098:                 energy=energy+PREFACTOR*C2T*SBMDHV(kappa,r1)1171:                 etemp=etemp+PREFACTOR*C2T*SBMDHV(kappa,r1)
1099: 1172: 
1100:                 f_over_r=PREFACTOR*C2T*SBMDHVP(kappa,r1)1173:                 f_over_r=PREFACTOR*C2T*SBMDHVP(kappa,r1)
1101:                 if(r2 .gt. DHswitch2)then1174:                 if(r2 .gt. DHswitch2)then
1102:                         RD1=r1-DHswitch1175:                 RD1=r1-DHswitch
1103:                         f_over_r=(f_over_r+C2T*(2*COEF2*RD1+3*COEF3*RD1**2))1176:                         f_over_r=(f_over_r+C2T*(2*COEF2*RD1+3*COEF3*RD1**2))
1104:                         energy=energy+COEF2*RD1**2+COEF3*RD1**31177:                         etemp=etemp+COEF2*RD1**2+COEF3*RD1**3
1105:                 elseif(r2 .lt. DHswitch2)then1178:                 elseif(r2 .lt. DHswitch2)then
1106:                 ! normal DH term1179:                 ! normal DH term
1107:                 else1180:                 else
1108:                 ! things should have fallen in one of the previous two...1181:                 ! things should have fallen in one of the previous two...
1109:                 write(*,*) 'something went wrong with DH switching function'1182:                 write(*,*) 'something went wrong with DH switching function'
1110: !                call abort1183: !                call abort
1111:                 endif1184:                 endif
1112:                 f_over_r=f_over_r*sqrt(rm2)1185:                 f_over_r=f_over_r*sqrt(rm2)
1113:               grad(C1*3-2) = grad(C1*3-2) + f_over_r * dx1186:               gradOMP(C1*3-2) = gradOMP(C1*3-2) + f_over_r * dx
1114:               grad(C1*3-1) = grad(C1*3-1) + f_over_r * dy1187:               gradOMP(C1*3-1) = gradOMP(C1*3-1) + f_over_r * dy
1115:               grad(C1*3)   = grad(C1*3)   + f_over_r * dz1188:               gradOMP(C1*3)   = gradOMP(C1*3)   + f_over_r * dz
1116:               grad(C2*3-2) = grad(C2*3-2) - f_over_r * dx1189: 
1117:               grad(C2*3-1) = grad(C2*3-1) - f_over_r * dy1190:               gradOMP(C2*3-2) =  gradOMP(C2*3-2) - f_over_r * dx
1118:               grad(C2*3)   = grad(C2*3)   - f_over_r * dz1191:               gradOMP(C2*3-1) =  gradOMP(C2*3-1) - f_over_r * dy
 1192:               gradOMP(C2*3)   =  gradOMP(C2*3)   - f_over_r * dz
1119: 1193: 
1120:             endif1194:             endif
1121: 1195: 
1122:            endif1196:            endif
1123:            enddo1197:            enddo
1124:  1198:  
1125:           endif1199:           endif
1126: 1200: 
1127:             enddo1201:             enddo
1128:            enddo1202:            enddo
1129:           enddo1203:           enddo
1130: 1204: 
1131:           enddo1205:           enddo
1132: !$OMP END PARALLEL1206: !$OMP END PARALLEL
 1207:                 energy = energy  + etemp
 1208:         do i=2,SBMCHARGEON(1)
 1209:            C1=SBMCHARGEON(i)
 1210:               grad(C1*3-2) = grad(C1*3-2)+ gradOMP(C1*3-2) 
 1211:               grad(C1*3-1) = grad(C1*3-1)+ gradOMP(C1*3-1) 
 1212:               grad(C1*3)   = grad(C1*3)  + gradOMP(C1*3)   
 1213:         enddo
 1214: 
1133:       END1215:       END
1134: 1216: 
1135: !^^^^^^^^^^^^^^^^^^^^^^^^^^^End of SBMDHELEC^^^^^^^^^^^^^^^^^^^^^^^^^1217: !^^^^^^^^^^^^^^^^^^^^^^^^^^^End of SBMDHELEC^^^^^^^^^^^^^^^^^^^^^^^^^
1136: 1218: 
1137: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<1219: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1138: !* SBMPR computes the forces due to position restraints               *1220: !* SBMPR computes the forces due to position restraints               *
1139: !**********************************************************************1221: !**********************************************************************
1140: 1222: 
1141:       subroutine SBMPR(x,y,z,grad, energy, natoms, 1223:       subroutine SBMPR(x,y,z,grad, energy, natoms, 
1142:      Q SBMPRN,SBMPRI,SBMPRK,SBMPRX)1224:      Q SBMPRN,SBMPRI,SBMPRK,SBMPRX)
1143: 1225: 
1144:       USE KEY1226:       USE KEY
1145:       implicit NONE1227:       implicit NONE
1146:       integer I, J, NATOMS1228:       integer I, J, NATOMS
1147:       INTEGER SBMPRN, SBMPRI(NATOMS)1229:       INTEGER SBMPRN, SBMPRI(NATOMS)
1148:       DOUBLE PRECISION SBMPRK(NATOMS,6),SBMPRX(NATOMS,3)1230:       DOUBLE PRECISION SBMPRK(NATOMS,6),SBMPRX(NATOMS,3)
1149:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS),1231:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS),
1150:      Q grad(3*NATOMS), energy1232:      Q grad(3*NATOMS), energy,GRADOMP(3*NATOMS), ETEMP
1151:       DOUBLE PRECISION DX,DY,DZ,K1,K2,K3,K4,K5,K61233:       DOUBLE PRECISION DX,DY,DZ,K1,K2,K3,K4,K5,K6
1152: 1234: 
1153: 1235: 
1154:       if(SBMPRN .eq. 0) RETURN1236:       if(SBMPRN .eq. 0) RETURN
1155: 1237: 
1156: !$OMP PARALLEL PRIVATE(I,J,DX,DY,DZ,K1,K2,K3,K4,K5,K6)REDUCTION(+:energy,grad)1238:       do I=1,NATOMS
 1239:          gradOMP(I*3)=0
 1240:          gradOMP(I*3-1)=0
 1241:          gradOMP(I*3-2)=0
 1242:       enddo
 1243: 
 1244: !$OMP PARALLEL PRIVATE(I,J,DX,DY,DZ,K1,K2,K3,K4,K5,K6)REDUCTION(+:etemp,gradOMP)
1157: !$OMP DO1245: !$OMP DO
1158:               do I=1,SBMPRN1246:               do I=1,SBMPRN
1159:            J=SBMPRI(I)1247:            J=SBMPRI(I)
1160:            DX=X(J)-SBMPRX(I,1)1248:            DX=X(J)-SBMPRX(I,1)
1161:            DY=Y(J)-SBMPRX(I,2)1249:            DY=Y(J)-SBMPRX(I,2)
1162:            DZ=Z(J)-SBMPRX(I,3)1250:            DZ=Z(J)-SBMPRX(I,3)
1163:            K1=SBMPRK(I,1)1251:            K1=SBMPRK(I,1)
1164:            K2=SBMPRK(I,2)1252:            K2=SBMPRK(I,2)
1165:            K3=SBMPRK(I,3)1253:            K3=SBMPRK(I,3)
1166:            K4=SBMPRK(I,4)1254:            K4=SBMPRK(I,4)
1167:            K5=SBMPRK(I,5)1255:            K5=SBMPRK(I,5)
1168:            K6=SBMPRK(I,6)1256:            K6=SBMPRK(I,6)
1169: 1257: 
1170: 1258: 
1171:            ENERGY = ENERGY +0.5*(K1*DX**2+K2*DY**2+K3*DZ**2) +1259:            ETEMP = ETEMP +0.5*(K1*DX**2+K2*DY**2+K3*DZ**2) +
1172:      Q     K4*DX*DY+K5*DX*DZ+K6*DY*DZ1260:      Q     K4*DX*DY+K5*DX*DZ+K6*DY*DZ
1173: 1261: 
1174:            grad(J*3-2) = grad(J*3-2) + K1*DX + K4*DY + K5*DZ 1262:            gradOMP(J*3-2) = gradOMP(J*3-2) + K1*DX + K4*DY + K5*DZ 
1175:            grad(J*3-1) = grad(J*3-1) + K2*DY + K4*DX + K6*DZ1263:            gradOMP(J*3-1) = gradOMP(J*3-1) + K2*DY + K4*DX + K6*DZ
1176:            grad(J*3)   = grad(J*3)   + K3*DZ + K5*DX + K6*DY1264:            gradOMP(J*3)   = gradOMP(J*3)   + K3*DZ + K5*DX + K6*DY
1177:         enddo1265:         enddo
1178: !$OMP END DO1266: !$OMP END DO
1179: !$OMP END PARALLEL1267: !$OMP END PARALLEL
 1268:         ENERGY=ENERGY+ETEMP
 1269:         do I=1,NATOMS
 1270:               grad(I*3-2) = grad(I*3-2)+ gradOMP(I*3-2) 
 1271:               grad(I*3-1) = grad(I*3-1)+ gradOMP(I*3-1) 
 1272:               grad(I*3)   = grad(I*3)  + gradOMP(I*3)   
 1273:         enddo
1180: 1274: 
1181:       END1275:       END
1182: 1276: 
1183: !^^^^^^^^^^^^^^^^^^^^^^^^^^^End of SBMPR^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^1277: !^^^^^^^^^^^^^^^^^^^^^^^^^^^End of SBMPR^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1184: 1278: 


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0