hdiff output

r30661/SBM.f 2016-07-06 15:35:44.619347138 +0100 r30660/SBM.f 2016-07-06 15:35:44.975351950 +0100
155:         read(30,*) ANtemp155:         read(30,*) ANtemp
156:         write(*,*) ANtemp, 'atoms'156:         write(*,*) ANtemp, 'atoms'
157:         NUMCHARGES=1157:         NUMCHARGES=1
158:         do i=1, ANtemp158:         do i=1, ANtemp
159:           read(30,*) TMPINT,TMPCHAR,TMPINT,TMPCHAR,TMPCHAR,SBMCHARGE(i),ATMASS(I)159:           read(30,*) TMPINT,TMPCHAR,TMPINT,TMPCHAR,TMPCHAR,SBMCHARGE(i),ATMASS(I)
160:           if(SBMCHARGE(i) .ne. 0)then160:           if(SBMCHARGE(i) .ne. 0)then
161:              NUMCHARGES=NUMCHARGES+1161:              NUMCHARGES=NUMCHARGES+1
162:              SBMCHARGEON(NUMCHARGES)=i162:              SBMCHARGEON(NUMCHARGES)=i
163:           endif 163:           endif 
164:         end do164:         end do
165:         SBMCHARGEON(1)=NUMCHARGES165:         SBMCHARGEON(1)=NUMCHARGES-1
166:         166:         
167:         read(30,*) 167:         read(30,*) 
168:         read(30,*) NC168:         read(30,*) NC
169:         write(*,*) NC, 'contacts'        169:         write(*,*) NC, 'contacts'        
170:         read(30,*) CONTACTTYPE170:         read(30,*) CONTACTTYPE
171:           if(NC .gt. MaxCon)then171:           if(NC .gt. MaxCon)then
172:              write(*,*) 'too many contacts'172:              write(*,*) 'too many contacts'
173:              STOP173:              STOP
174:           endif174:           endif
175: 175: 
574:             z10 = 1.0/FXI574:             z10 = 1.0/FXI
575:             z20 = 1.0/FYI575:             z20 = 1.0/FYI
576:             Z12 = Z10*Z20576:             Z12 = Z10*Z20
577:             Z1 = Z10577:             Z1 = Z10
578:             Z2 = Z20578:             Z2 = Z20
579:             CT0 = MIN(1.0,CT*Z12)579:             CT0 = MIN(1.0,CT*Z12)
580:             CT1 = MAX(-1.0,CT0)580:             CT1 = MAX(-1.0,CT0)
581:             S = XKJ*(DZ*GY-DY*GZ)+YKJ*(DX*GZ-DZ*GX)+ZKJ*(DY*GX-DX*GY)581:             S = XKJ*(DZ*GY-DY*GZ)+YKJ*(DX*GZ-DZ*GX)+ZKJ*(DY*GX-DX*GY)
582:             AP0 = ACOS(CT1)582:             AP0 = ACOS(CT1)
583:             AP1 = PI-SIGN(AP0,S)583:             AP1 = PI-SIGN(AP0,S)
 584: 
 585: 
584:             CT = AP1586:             CT = AP1
585:             CPHI = COS(AP1)587:             CPHI = COS(AP1)
586:             SPHI = SIN(AP1)588:             SPHI = SIN(AP1)
 589: 
 590: 
 591: !            CT0 = CT(II)
587: ! Here is the energy part592: ! Here is the energy part
588:           A1=CT-PHISBM(II)593:           A1=CT-PHISBM(II)
589:           A3=A1*3594:           A3=A1*3
590: 595: 
591:         if(PHITYPE(II) .eq. 1)then596:         if(PHITYPE(II) .eq. 1)then
592:           etemp =  etemp + PK(II)*(3.0/2.0-cos(A1)-0.5*cos(A3))597:           etemp =  etemp + PK(II)*(3.0/2.0-cos(A1)-0.5*cos(A3))
593: ! dE/dPHI598: ! dE/dPHI
594:             DF=PK(II)*(sin(A1)+1.5*sin(A3))599:             DF=PK(II)*(sin(A1)+1.5*sin(A3))
595:         elseif(PHITYPE(II) .eq. 2)then600:         elseif(PHITYPE(II) .eq. 2)then
596:           if(A1 .gt. PI)then601:           if(A1 .gt. PI)then
607: 612: 
608:         write(*,*) 'unrecognized dihedral type', PHITYPE(II)613:         write(*,*) 'unrecognized dihedral type', PHITYPE(II)
609:         STOP614:         STOP
610:         endif615:         endif
611: 616: 
612: ! insert the new 617: ! insert the new 
613: 618: 
614: ! now, do dPhi/dX619: ! now, do dPhi/dX
615: 620: 
616: ! |G|/|A|**2 621: ! |G|/|A|**2 
617:           TT1 = Z1*Z1*RKJ*DF622:             TT1 = Z1*Z1*RKJ*DF
618: ! FG/(A**2*|G|)623: ! FG/(A**2*|G|)
619:           TT2 = FGoverG*Z1*Z1*DF624:             TT2 = FGoverG*Z1*Z1*DF
620: ! HG/(B**2*|G|)625: ! HG/(B**2*|G|)
621:           TT3 = HGoverG*Z2*Z2*DF626:             TT3 = HGoverG*Z2*Z2*DF
622: ! |G|/|B|**2 627: ! |G|/|B|**2 
623:           TT4 = Z2*Z2*RKJ*DF628:             TT4 = Z2*Z2*RKJ*DF
624: ! note: negatives are flipped from paper because A=-DX629: ! note: negatives are flipped from paper because A=-DX
625:           TT1X=TT1*DX630:         TT1X=TT1*DX
626:           TT1Y=TT1*DY631:         TT1Y=TT1*DY
627:           TT1Z=TT1*DZ632:         TT1Z=TT1*DZ
628:           TT2X=TT2*DX633: 
629:           TT2Y=TT2*DY634:         TT2X=TT2*DX
630:           TT2Z=TT2*DZ635:         TT2Y=TT2*DY
631:           TT3X=TT3*GX636:         TT2Z=TT2*DZ
632:           TT3Y=TT3*GY637: 
633:           TT3Z=TT3*GZ638: 
634:           TT4X=TT4*GX639:         TT3X=TT3*GX
635:           TT4Y=TT4*GY640:         TT3Y=TT3*GY
636:           TT4Z=TT4*GZ641:         TT3Z=TT3*GZ
637:           gradOMP(I3*3-2) =  gradOMP(I3*3-2)  + TT1X  642: 
638:           gradOMP(I3*3-1) =  gradOMP(I3*3-1)  + TT1Y643: 
639:           gradOMP(I3*3)   =  gradOMP(I3*3)    + TT1Z644:         TT4X=TT4*GX
640:           gradOMP(J3*3-2) =  gradOMP(J3*3-2)  - TT1X - TT2X - TT3X645:         TT4Y=TT4*GY
641:           gradOMP(J3*3-1) =  gradOMP(J3*3-1)  - TT1Y - TT2Y - TT3Y646:         TT4Z=TT4*GZ
642:           gradOMP(J3*3)   =  gradOMP(J3*3)    - TT1Z - TT2Z - TT3Z647: 
643:           gradOMP(K3*3-2) =  gradOMP(K3*3-2)  + TT2X + TT3X - TT4X648:             gradOMP(I3*3-2) =  gradOMP(I3*3-2)  + TT1X  
644:           gradOMP(K3*3-1) =  gradOMP(K3*3-1)  + TT2Y + TT3Y - TT4Y649:             gradOMP(I3*3-1) =  gradOMP(I3*3-1)  + TT1Y
645:           gradOMP(K3*3)   =  gradOMP(K3*3)    + TT2Z + TT3Z - TT4Z650:             gradOMP(I3*3)   =  gradOMP(I3*3)    + TT1Z
646:           gradOMP(L3*3-2) =  gradOMP(L3*3-2)  + TT4X651:             gradOMP(J3*3-2) =  gradOMP(J3*3-2)  - TT1X - TT2X - TT3X
647:           gradOMP(L3*3-1) =  gradOMP(L3*3-1)  + TT4Y652:             gradOMP(J3*3-1) =  gradOMP(J3*3-1)  - TT1Y - TT2Y - TT3Y
648:           gradOMP(L3*3)   =  gradOMP(L3*3)    + TT4Z653:             gradOMP(J3*3)   =  gradOMP(J3*3)    - TT1Z - TT2Z - TT3Z
 654:             gradOMP(K3*3-2) =  gradOMP(K3*3-2)  + TT2X + TT3X - TT4X
 655:             gradOMP(K3*3-1) =  gradOMP(K3*3-1)  + TT2Y + TT3Y - TT4Y
 656:             gradOMP(K3*3)   =  gradOMP(K3*3)    + TT2Z + TT3Z - TT4Z
 657:             gradOMP(L3*3-2) =  gradOMP(L3*3-2)  + TT4X
 658:             gradOMP(L3*3-1) =  gradOMP(L3*3-1)  + TT4Y
 659:             gradOMP(L3*3)   =  gradOMP(L3*3)    + TT4Z
649:         END DO660:         END DO
650: !$OMP END DO 661: !$OMP END DO 
651: !$OMP END PARALLEL662: !$OMP END PARALLEL
652: 663: 
653:         ENERGY=ENERGY+etemp664:         ENERGY=ENERGY+etemp
654:         do C1=1,NATOMS665:         do C1=1,NATOMS
655:               grad(C1*3-2) = grad(C1*3-2)+ gradOMP(C1*3-2) 666:               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) 667:               grad(C1*3-1) = grad(C1*3-1)+ gradOMP(C1*3-1) 
657:               grad(C1*3)   = grad(C1*3)  + gradOMP(C1*3)   668:               grad(C1*3)   = grad(C1*3)  + gradOMP(C1*3)   
658:         enddo669:         enddo
666: 677: 
667: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<678: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
668: !* SBMCONTACTS: computes the force on all atoms due to contacts via a   *679: !* SBMCONTACTS: computes the force on all atoms due to contacts via a   *
669: !* 10-12 or 6-12 potential                                              *680: !* 10-12 or 6-12 potential                                              *
670: !***********************************************************************681: !***********************************************************************
671: 682: 
672:       subroutine SBMcontacts(x,y,z,grad,energy,683:       subroutine SBMcontacts(x,y,z,grad,energy,
673:      Q NATOMS,IC,JC,Sigma,EpsC,NC,CONTACTTYPE)684:      Q NATOMS,IC,JC,Sigma,EpsC,NC,CONTACTTYPE)
674:       USE KEY685:       USE KEY
675:       implicit NONE686:       implicit NONE
676:       integer I, NATOMS,NC687:       integer I, N, J,NATOMS,NC
677: 688: 
678:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS), 689:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS) 
679:      Q grad(3*NATOMS), energy,etemp, GRADOMP(3*NATOMS)690:      Q , grad(3*NATOMS), energy
680:       DOUBLE PRECISION dx,dy,dz691:       DOUBLE PRECISION dx,dy,dz
681: 692: 
682:       integer C1, C2 693:       integer C1, C2, ConfID, Q, SC1, SC2, Cf1
683:       DOUBLE PRECISION  r2, rm2, rm10, rm6,f_over_r,r1 694:       DOUBLE PRECISION  r2, rm2, rm10, rm6,f_over_r, dsig, deps, 
 695:      Q s1, s2, ep1, ep2, r1, rc,r, summm
684: 696: 
685:         DOUBLE PRECISION Sigma(NC), EpsC(NC)697:         DOUBLE PRECISION Sigma(NC), EpsC(NC)
686:         INTEGER IC(NC), JC(NC),CONTACTTYPE698:         INTEGER IC(NC), JC(NC),CONTACTTYPE
687: 699: 
688: ! type 1 is 6-12 
689: ! type 2 is 12-12 
690: ! implement type 3 as gaussian700: ! implement type 3 as gaussian
691: ! implement type 4 as dual-basin gaussian701: ! implement type 4 as dual-basin gaussian
692: ! question: Should type be part of the interaction, or should they be organized702: ! question: Should type be part of the interaction, or should they be organized
693: ! into separate lists?703: ! 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: 704: 
700: 705: 
701:         ETEMP=0 
702: ! type 1 is 6-12 interaction706: ! type 1 is 6-12 interaction
703:         if(CONTACTTYPE .eq. 1)then707:         if(CONTACTTYPE .eq. 1)then
704: !$OMP PARALLEL PRIVATE(I,C1,C2,DX,DY,DZ,R2,RM2,RM6,F_OVER_R)REDUCTION(+:etemp,gradOMP)708: 
705: !$OMP DO709:                 do i=1, NC
706:           do i=1, NC710:                
707:             C1 = IC(i)711:                         C1 = IC(i)
708:             C2 = JC(i)712:                         C2 = JC(i)
709:             dx = X(C1) - X(C2)713:                         dx = X(C1) - X(C2)
710:             dy = Y(C1) - Y(C2)714:                         dy = Y(C1) - Y(C2)
711:             dz = Z(C1) - Z(C2)715:                         dz = Z(C1) - Z(C2)
712:             r2 = dx**2 + dy**2 + dz**2716:                 
713:             rm2 = 1.0/r2717:                         r2 = dx**2 + dy**2 + dz**2
714:             rm2 = rm2*sigma(i)718:                         rm2 = 1.0/r2
715:             rm6 = rm2**3719:                         rm2 = rm2*sigma(i)
716:             ETEMP = ETEMP + epsC(i)*rm6*(rm6-2.0)720:                         rm6 = rm2**3
717:             f_over_r = -epsC(i)*12.0*rm6*(rm6-1.0)/r2721:                 
718:             gradOMP(3*C1-2) = gradOMP(3*C1-2) + f_over_r * dx722:                         energy = energy + epsC(i)*rm6*(rm6-2.0)
719:             gradOMP(3*C1-1) = gradOMP(3*C1-1) + f_over_r * dy723:                 
720:             gradOMP(3*C1)   = gradOMP(3*C1)   + f_over_r * dz724:                         f_over_r = -epsC(i)*12.0*rm6*(rm6-1.0)/r2
721:             gradOMP(3*C2-2) = gradOMP(3*C2-2) - f_over_r * dx725:                 
722:             gradOMP(3*C2-1) = gradOMP(3*C2-1) - f_over_r * dy726:                         grad(3*C1-2) = grad(3*C1-2) + f_over_r * dx
723:             gradOMP(3*C2)   = gradOMP(3*C2)   - f_over_r * dz727:                         grad(3*C1-1) = grad(3*C1-1) + f_over_r * dy
724:           enddo728:                         grad(3*C1)   = grad(3*C1)   + f_over_r * dz
725: !$OMP END DO729:                 
726: !$OMP END PARALLEL730:                         grad(3*C2-2) =  grad(3*C2-2) - f_over_r * dx
 731:                         grad(3*C2-1) =  grad(3*C2-1) - f_over_r * dy
 732:                         grad(3*C2)   =  grad(3*C2)   - f_over_r * dz
 733:                 enddo
 734: 
727: ! type 2 is 10-12 interaction735: ! type 2 is 10-12 interaction
728:         elseif(CONTACTTYPE .eq. 2)then736:         elseif(CONTACTTYPE .eq. 2)then
729: !$OMP PARALLEL PRIVATE(I,C1,C2,DX,DY,DZ,R2,RM2,RM10,F_OVER_R)REDUCTION(+:etemp,gradOMP)737:         do i=1, NC
730: !$OMP DO738:        
731:           do i=1, NC739:                 C1 = IC(i)
732:              C1 = IC(i)740:                 C2 = JC(i)
733:              C2 = JC(i)741:                 dx = X(C1) - X(C2)
734:              dx = X(C1) - X(C2)742:                 dy = Y(C1) - Y(C2)
735:              dy = Y(C1) - Y(C2)743:                 dz = Z(C1) - Z(C2)
736:              dz = Z(C1) - Z(C2)744:         
737:              r2 = dx**2 + dy**2 + dz**2745:                 r2 = dx**2 + dy**2 + dz**2
738:              rm2 = 1.0/r2746:         
739:              rm2 = rm2*sigma(i)747:                 rm2 = 1.0/r2
740:              rm10 = rm2**5748:                 rm2 = rm2*sigma(i)
741:              ETEMP = ETEMP + epsC(i)*rm10*(5*rm2-6.0)749:                 rm10 = rm2**5
742:              f_over_r = -epsc(i)*60.0*rm10*(rm2-1.0)/r2750:                 
743:              gradOMP(3*C1-2) = gradOMP(3*C1-2) + f_over_r * dx751:                 energy = energy + epsC(i)*rm10*(5*rm2-6.0)
744:              gradOMP(3*C1-1) = gradOMP(3*C1-1) + f_over_r * dy752:                 f_over_r = -epsc(i)*60.0*rm10*(rm2-1.0)/r2
745:              gradOMP(3*C1)   = gradOMP(3*C1)   + f_over_r * dz753:                 
746:              gradOMP(3*C2-2) = gradOMP(3*C2-2) - f_over_r * dx754:                 grad(3*C1-2) = grad(3*C1-2) + f_over_r * dx
747:              gradOMP(3*C2-1) = gradOMP(3*C2-1) - f_over_r * dy755:                 grad(3*C1-1) = grad(3*C1-1) + f_over_r * dy
748:              gradOMP(3*C2)   = gradOMP(3*C2)   - f_over_r * dz756:                 grad(3*C1)   = grad(3*C1)   + f_over_r * dz
749:            enddo757:                 
750: !$OMP END DO758:                 grad(3*C2-2) =  grad(3*C2-2) - f_over_r * dx
751: !$OMP END PARALLEL759:                 grad(3*C2-1) =  grad(3*C2-1) - f_over_r * dy
 760:                 grad(3*C2)   =  grad(3*C2)   - f_over_r * dz
 761:               enddo
752: 762: 
753:         else763:         else
754:         write(*,*) CONTACTTYPE, 'is not a valid contact selection'764:         write(*,*) CONTACTTYPE, 'is not a valid contact selection'
755:         stop765:         stop
756:         endif766:         endif
757: 767: 
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:  
766: 768: 
767:       end769:       end
768: 770: 
769: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^end of SBMContacts^^^^^^^^^^^^^^^^^^^^^^^^^^^771: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^end of SBMContacts^^^^^^^^^^^^^^^^^^^^^^^^^^^
770: 772: 
771: 773: 
772: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<774: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
773: !* SBMNonContacts computes the forces due to non native contacts       *775: !* SBMNonContacts computes the forces due to non native contacts       *
774: !**********************************************************************776: !**********************************************************************
775: 777: 
1011:      Q grad(3*NATOMS), gradOMP(3*NATOMS), energy,STT,SA,SB,SC,1013:      Q grad(3*NATOMS), gradOMP(3*NATOMS), energy,STT,SA,SB,SC,
1012:      Q SBMDHVP, SBMDHV1014:      Q SBMDHVP, SBMDHV
1013:       integer C1, C2, ii,jj,kk,k,l,iii,jjj1015:       integer C1, C2, ii,jj,kk,k,l,iii,jjj
1014:       DOUBLE PRECISION  r1,r2, rm2, rm14, f_over_r, NCswitch,NCcut 1016:       DOUBLE PRECISION  r1,r2, rm2, rm14, f_over_r, NCswitch,NCcut 
1015:       DOUBLE PRECISION A,B,C, D, COEF2, COEF31017:       DOUBLE PRECISION A,B,C, D, COEF2, COEF3
1016:         INTEGER NTHREADS,  OMP_GET_NUM_THREADS,1018:         INTEGER NTHREADS,  OMP_GET_NUM_THREADS,
1017:      Q  OMP_GET_THREAD_NUM,OMP_NUM_THREADS1019:      Q  OMP_GET_THREAD_NUM,OMP_NUM_THREADS
1018: 1020: 
1019:         INTEGER NNCinc(NATOMS,NATOMS)1021:         INTEGER NNCinc(NATOMS,NATOMS)
1020:         DOUBLE PRECISION dx,dy,dz1022:         DOUBLE PRECISION dx,dy,dz
1021:         integer tempN, alpha,SBMCHARGEON(NATOMS+1)1023:         integer tempN, alpha,SBMCHARGEON(NATOMS)
1022:         double precision diff,Vfunc,Ffunc,SBMCHARGE(NATOMS)1024:         double precision diff,Vfunc,Ffunc,SBMCHARGE(NATOMS)
1023:         double precision PREFACTOR,KAPPA,DHswitch,DHswitch2,DHcut,DHcut21025:         double precision PREFACTOR,KAPPA,DHswitch,DHswitch2,DHcut,DHcut2
1024:         ! Ngrid is the number of atoms in that grid point1026:         ! Ngrid is the number of atoms in that grid point
1025:         ! grid is the array of atoms in each grid1027:         ! grid is the array of atoms in each grid
1026:         integer Ngrid,grid,maxgrid,maxpergrid1028:         integer Ngrid,grid,maxgrid,maxpergrid
1027:         ! number of atoms per grid, max1029:         ! number of atoms per grid, max
1028:         parameter (maxpergrid=50)1030:         parameter (maxpergrid=50)
1029:         ! dimensions of grid1031:         ! dimensions of grid
1030:         parameter (maxgrid=15)1032:         parameter (maxgrid=15)
1031:         dimension Ngrid(maxgrid,maxgrid,maxgrid),1033:         dimension Ngrid(maxgrid,maxgrid,maxgrid),
1032:      Q  grid(maxgrid,maxgrid,maxgrid,maxpergrid)1034:      Q  grid(maxgrid,maxgrid,maxgrid,maxpergrid)
1033:         integer MaxGridX,MaxGridY,MaxGridZ1035:         integer MaxGridX,MaxGridY,MaxGridZ
1034:         double precision gridsize,RD11036:         double precision gridsize,RD1
1035:         double precision minX,minY,minZ,maxX,maxY,maxZ1037:         double precision minX,minY,minZ,maxX,maxY,maxZ
1036:         integer Xgrid,Ygrid,Zgrid,TID1038:         integer Xgrid,Ygrid,Zgrid,TID
1037:         double precision Etemp,C2T1039:         double precision Etemp,C2T
1038:  
1039:  
1040:       if(SBMCHARGEON(1) .eq. 1) RETURN 
1041:  
1042:         diff=DHcut-DHswitch1040:         diff=DHcut-DHswitch
1043:         alpha=121041:         alpha=12
1044: 1042: 
1045: 1043: 
1046:         GRIDSIZE=DHCUT*1.011044:         GRIDSIZE=DHCUT*1.01
1047:         DHcut2=DHcut**21045:         DHcut2=DHcut**2
1048:         DHswitch2=DHswitch**21046:         DHswitch2=DHswitch**2
1049: 1047: 
1050:         A=-PREFACTOR*SBMDHV(kappa,DHcut)1048:         A=-PREFACTOR*SBMDHV(kappa,DHcut)
1051:         B=-PREFACTOR*SBMDHVP(kappa,DHcut)1049:         B=-PREFACTOR*SBMDHVP(kappa,DHcut)
1222: 1220: 
1223:       subroutine SBMPR(x,y,z,grad, energy, natoms, 1221:       subroutine SBMPR(x,y,z,grad, energy, natoms, 
1224:      Q SBMPRN,SBMPRI,SBMPRK,SBMPRX)1222:      Q SBMPRN,SBMPRI,SBMPRK,SBMPRX)
1225: 1223: 
1226:       USE KEY1224:       USE KEY
1227:       implicit NONE1225:       implicit NONE
1228:       integer I, J, NATOMS1226:       integer I, J, NATOMS
1229:       INTEGER SBMPRN, SBMPRI(NATOMS)1227:       INTEGER SBMPRN, SBMPRI(NATOMS)
1230:       DOUBLE PRECISION SBMPRK(NATOMS,6),SBMPRX(NATOMS,3)1228:       DOUBLE PRECISION SBMPRK(NATOMS,6),SBMPRX(NATOMS,3)
1231:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS),1229:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS),
1232:      Q grad(3*NATOMS), energy,GRADOMP(3*NATOMS), ETEMP1230:      Q grad(3*NATOMS), energy
1233:       DOUBLE PRECISION DX,DY,DZ,K1,K2,K3,K4,K5,K61231:       DOUBLE PRECISION DX,DY,DZ,K1,K2,K3,K4,K5,K6
1234:  
1235:  
1236:       if(SBMPRN .eq. 0) RETURN 
1237:  
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) 
1245: !$OMP DO 
1246:               do I=1,SBMPRN1232:               do I=1,SBMPRN
1247:            J=SBMPRI(I)1233:            J=SBMPRI(I)
1248:            DX=X(J)-SBMPRX(I,1)1234:            DX=X(J)-SBMPRX(I,1)
1249:            DY=Y(J)-SBMPRX(I,2)1235:            DY=Y(J)-SBMPRX(I,2)
1250:            DZ=Z(J)-SBMPRX(I,3)1236:            DZ=Z(J)-SBMPRX(I,3)
1251:            K1=SBMPRK(I,1)1237:            K1=SBMPRK(I,1)
1252:            K2=SBMPRK(I,2)1238:            K2=SBMPRK(I,2)
1253:            K3=SBMPRK(I,3)1239:            K3=SBMPRK(I,3)
1254:            K4=SBMPRK(I,4)1240:            K4=SBMPRK(I,4)
1255:            K5=SBMPRK(I,5)1241:            K5=SBMPRK(I,5)
1256:            K6=SBMPRK(I,6)1242:            K6=SBMPRK(I,6)
1257: 1243: 
1258: 1244: 
1259:            ETEMP = ETEMP +0.5*(K1*DX**2+K2*DY**2+K3*DZ**2) +1245:            ENERGY = ENERGY + 
1260:      Q     K4*DX*DY+K5*DX*DZ+K6*DY*DZ1246:      Q     0.5*K1*DX**2 + 
1261: 1247:      Q     0.5*K2*DY**2 +
1262:            gradOMP(J*3-2) = gradOMP(J*3-2) + K1*DX + K4*DY + K5*DZ 1248:      Q     0.5*K3*DZ**2 +
1263:            gradOMP(J*3-1) = gradOMP(J*3-1) + K2*DY + K4*DX + K6*DZ1249:      Q     K4*DX*DY +
1264:            gradOMP(J*3)   = gradOMP(J*3)   + K3*DZ + K5*DX + K6*DY1250:      Q     K5*DX*DZ +
1265:         enddo1251:      Q     K6*DY*DZ
1266: !$OMP END DO1252: 
1267: !$OMP END PARALLEL1253:            grad(J*3-2) = grad(J*3-2) + K1*DX + K4*DY + K5*DZ 
1268:         ENERGY=ENERGY+ETEMP1254:            grad(J*3-1) = grad(J*3-1) + K2*DY + K4*DX + K6*DZ
1269:         do I=1,NATOMS1255:            grad(J*3)   = grad(J*3)   + K3*DZ + K5*DX + K6*DY
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:         enddo1256:         enddo
1274: 1257: 
1275:       END1258:       END
1276: 1259: 
1277: !^^^^^^^^^^^^^^^^^^^^^^^^^^^End of SBMPR^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^1260: !^^^^^^^^^^^^^^^^^^^^^^^^^^^End of SBMPR^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1278: 1261: 


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0