criticall bug fix in energy
[unres4.git] / source / unres / energy.F90
index 4e043fe..3870bb0 100644 (file)
 !          print *,"Processor",myrank," BROADCAST iorder"
 ! FG master sets up the WEIGHTS_ array which will be broadcast to the 
 ! FG slaves as WEIGHTS array.
-         ! weights_(1)=wsc
+          weights_(1)=wsc
           weights_(2)=wscp
           weights_(3)=welec
           weights_(4)=wcorr
              .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.d0 &
              .or. wcorr6.gt.0.0d0 .or. wturn6.gt.0.0d0 ) then
 #endif
-!            write(iout,*),"just befor eelec call"
+!            print *,"just befor eelec call"
             call eelec(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
-!         write (iout,*) "ELEC calc"
+!            print *, "ELEC calc"
          else
             ees=0.0d0
             evdw1=0.0d0
        eespp=0.0d0
       endif
 !      write(iout,*) ecorr_nucl,"ecorr_nucl",nres_molec(2)
+!      print *,"before ecatcat"
       if (nfgtasks.gt.1) then
       if (fg_rank.eq.0) then
       call ecatcat(ecationcation)
 ! to calculate the el-loc multibody terms of various order.
 !
 !AL el      mu=0.0d0
+   
 #ifdef PARMAT
       do i=ivec_start+2,ivec_end+2
 #else
       do i=3,nres+1
 #endif
         if (i.gt. nnt+2 .and. i.lt.nct+2) then
+          if (itype(i-2,1).eq.0) then 
+          iti = nloctyp
+          else
           iti = itype2loc(itype(i-2,1))
+          endif
         else
           iti=nloctyp
         endif
 !d        write (iout,*) 'mu2',mu2(:,i-2)
         if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.wcorr6.gt.0.0d0) &
         then  
-        call matmat2(CC(1,1,i-2),Ugder(1,1,i-2),CUgder(1,1,i-2))
+        call matmat2(CC(1,1,i-1),Ugder(1,1,i-2),CUgder(1,1,i-2))
         call matmat2(DD(1,1,i-2),Ugder(1,1,i-2),DUgder(1,1,i-2))
         call matmat2(Dtilde(1,1,i-2),Ug2der(1,1,i-2),DtUg2der(1,1,i-2))
-        call matvec2(Ctilde(1,1,i-2),obrot_der(1,i-2),Ctobrder(1,i-2))
+        call matvec2(Ctilde(1,1,i-1),obrot_der(1,i-2),Ctobrder(1,i-2))
         call matvec2(Dtilde(1,1,i-2),obrot2_der(1,i-2),Dtobr2der(1,i-2))
 ! Vectors and matrices dependent on a single virtual-bond dihedral.
         call matvec2(DD(1,1,i-2),b1tilde(1,iti1),auxvec(1))
         call matvec2(Ug2(1,1,i-2),auxvec(1),Ug2Db1t(1,i-2)) 
         call matvec2(Ug2der(1,1,i-2),auxvec(1),Ug2Db1tder(1,i-2)) 
-        call matvec2(CC(1,1,i-2),Ub2(1,i-2),CUgb2(1,i-2))
-        call matvec2(CC(1,1,i-2),Ub2der(1,i-2),CUgb2der(1,i-2))
+        call matvec2(CC(1,1,i-1),Ub2(1,i-2),CUgb2(1,i-2))
+        call matvec2(CC(1,1,i-1),Ub2der(1,i-2),CUgb2der(1,i-2))
         call matmat2(EUg(1,1,i-2),CC(1,1,iti1),EUgC(1,1,i-2))
         call matmat2(EUgder(1,1,i-2),CC(1,1,iti1),EUgCder(1,1,i-2))
         call matmat2(EUg(1,1,i-2),DD(1,1,iti1),EUgD(1,1,i-2))
           +a23*gmuij1(2)&
           +a32*gmuij1(3)&
           +a33*gmuij1(4))&
-         *fac_shield(i)*fac_shield(j)
+         *fac_shield(i)*fac_shield(j)&
+                    *sss_ele_cut
+
 !c         write(iout,*) "derivative over thatai"
 !c         write(iout,*) a22*gmuij1(1), a23*gmuij1(2) ,a32*gmuij1(3),
 !c     &   a33*gmuij1(4) 
           +a33*gmuij2(4)
          gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+&
            geel_loc_ij*wel_loc&
-         *fac_shield(i)*fac_shield(j)
+         *fac_shield(i)*fac_shield(j)&
+                    *sss_ele_cut
+
 
 !c  Derivative over j residue
          geel_loc_ji=a22*gmuji1(1)&
 
         gloc(nphi+j,icg)=gloc(nphi+j,icg)+&
            geel_loc_ji*wel_loc&
-         *fac_shield(i)*fac_shield(j)
+         *fac_shield(i)*fac_shield(j)&
+                    *sss_ele_cut
+
 
          geel_loc_ji=&
           +a22*gmuji2(1)&
 !c     &   a33*gmuji2(4)
          gloc(nphi+j-1,icg)=gloc(nphi+j-1,icg)+&
            geel_loc_ji*wel_loc&
-         *fac_shield(i)*fac_shield(j)
+         *fac_shield(i)*fac_shield(j)&
+                    *sss_ele_cut
 #endif
 
 !          write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
 #ifdef TIMING
       time01=MPI_Wtime()
 #endif
+!#define DEBUG
 #ifdef DEBUG
       write (iout,*) "sum_gradient gvdwc, gvdwx"
       do i=1,nres
 !        write (iout,'(i5,3f10.5,2x,f10.5)') 
 !     &  i,(gcorr4_turn(j,i),j=1,3),gel_loc_turn4(i)
 !      enddo
-      write (iout,*) "gvdwc gvdwc_scp gvdwc_scpp"
-      do i=1,nres
-        write (iout,'(i3,3f10.5,5x,3f10.5,5x,f10.5)') &
-         i,(gvdwc(j,i),j=1,3),(gvdwc_scp(j,i),j=1,3),&
-         (gvdwc_scpp(j,i),j=1,3)
-      enddo
-      write (iout,*) "gelc_long gvdwpp gel_loc_long"
-      do i=1,nres
-        write (iout,'(i3,3f10.5,5x,3f10.5,5x,f10.5)') &
-         i,(gelc_long(j,i),j=1,3),(gvdwpp(j,i),j=1,3),&
-         (gelc_loc_long(j,i),j=1,3)
-      enddo
+!      write (iout,*) "gvdwc gvdwc_scp gvdwc_scpp"
+!      do i=1,nres
+!        write (iout,'(i3,3f10.5,5x,3f10.5,5x,f10.5)') &
+!         i,(gvdwc(j,i),j=1,3),(gvdwc_scp(j,i),j=1,3),&
+!         (gvdwc_scpp(j,i),j=1,3)
+!      enddo
+!      write (iout,*) "gelc_long gvdwpp gel_loc_long"
+!      do i=1,nres
+!        write (iout,'(i3,3f10.5,5x,3f10.5,5x,f10.5)') &
+!         i,(gelc_long(j,i),j=1,3),(gvdwpp(j,i),j=1,3),&
+!         (gelc_loc_long(j,i),j=1,3)
+!      enddo
       call flush(iout)
 #endif
 #ifdef SPLITELE
 !      call intcartderiv
 !      call checkintcartgrad
       call zerograd
-      aincr=2.0D-5
+      aincr=1.0D-7
       write(iout,*) 'Calling CHECK_ECARTINT.',aincr
       nf=0
       icall=0
         rs,rthrp,rfourp,rsixp,reight,Irsixp,Ireight,Irtw,Irfourt,      &
         opt1,opt2,opt3,opt4,opt5,opt6,opt7,opt8,opt9,opt10,opt11,opt12,&
         opt13,opt14,opt15,opt16,opt17,opt18,opt19, &
-        Equad1,Equad2,dscmag,v1dpv2,dscmag3,constA,constB,Edip
+        Equad1,Equad2,dscmag,v1dpv2,dscmag3,constA,constB,Edip,&
+        ndiv,ndivi
         real(kind=8),dimension(3) ::dEvan1Cmcat,dEvan2Cmcat,dEeleccat,&
         gg,r,EtotalCat,dEtotalCm,dEtotalCalp,dEvan1Cm,dEvan2Cm, &
         dEtotalpep,dEtotalcat_num,dEddci,dEtotalcm_num,dEtotalcalp_num, &
         ecation_prot=0.0d0
 ! first lets calculate interaction with peptide groups
         if (nres_molec(5).eq.0) return
-         wconst=78
-        wdip =1.092777950857032D2
-        wdip=wdip/wconst
-        wmodquad=-2.174122713004870D4
-        wmodquad=wmodquad/wconst
-        wquad1 = 3.901232068562804D1
-        wquad1=wquad1/wconst
-        wquad2 = 3
-        wquad2=wquad2/wconst
-        wvan1 = 0.1
-        wvan2 = 6
         itmp=0
         do i=1,4
         itmp=itmp+nres_molec(i)
           if (zi.lt.0) zi=zi+boxzsize
 
          do j=itmp+1,itmp+nres_molec(5)
+!           print *,"WTF",itmp,j,i
+! all parameters were for Ca2+ to approximate single charge divide by two
+         ndiv=1.0
+         if ((itype(j,5).eq.1).or.(itype(j,5).eq.3)) ndiv=2.0
+         wconst=78*ndiv
+        wdip =1.092777950857032D2
+        wdip=wdip/wconst
+        wmodquad=-2.174122713004870D4
+        wmodquad=wmodquad/wconst
+        wquad1 = 3.901232068562804D1
+        wquad1=wquad1/wconst
+        wquad2 = 3
+        wquad2=wquad2/wconst
+        wvan1 = 0.1
+        wvan2 = 6
+!        itmp=0
+
            xj=c(1,j)
            yj=c(2,j)
            zj=c(3,j)
         E2 = -wquad1*Irthrp*wquad2+wvan1*(wvan2**12*Irtwelv-&
              2*wvan2**6*Irsistp)
         ecation_prot = ecation_prot+E1+E2
+!        print *,"ecatprot",i,j,ecation_prot,rcpm
         dE1dr = -2*costhet*wdip*Irthrp-& 
          (4*wmodquad*Irfiftp+3*wquad1*Irfourp)*sin2thet
         dE2dr = 3*wquad1*wquad2*Irfourp-     &
           enddo
            cm1mag=sqrt(cm1(1)**2+cm1(2)**2+cm1(3)**2)
          do j=itmp+1,itmp+nres_molec(5)
+         ndiv=1.0
+         if ((itype(j,5).eq.1).or.(itype(j,5).eq.3)) ndiv=2.0
+
            xj=c(1,j)
            yj=c(2,j)
            zj=c(3,j)
        endif
 !       enddo
 !       enddo
-         if(itype(i,1).eq.15.or.itype(i,1).eq.16) then
+! 15- Glu 16-Asp
+         if((itype(i,1).eq.15.or.itype(i,1).eq.16).or.&
+         ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.&
+         (itype(i,1).eq.25))) then
             if(itype(i,1).eq.16) then
             inum=1
             else
 
 !  The weights of the energy function calculated from
 !The quantum mechanical GAMESS simulations of calcium with ASP/GLU
-        wh2o=78
+          if ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.(itype(i,1).eq.25)) then
+            ndivi=0.5
+          else
+            ndivi=1.0
+          endif
+         ndiv=1.0
+         if ((itype(j,5).eq.1).or.(itype(j,5).eq.3)) ndiv=2.0
+
+        wh2o=78*ndivi*ndiv
         wc = vcatprm(1)
         wc=wc/wh2o
         wdip =vcatprm(2)
         v1dpv2 = v1(1)*v2(1)+v1(2)*v2(2)+v1(3)*v2(3)
 !  The weights of the energy function calculated from
 !The quantum mechanical GAMESS simulations of ASN/GLN with calcium
-        wh2o=78
+         ndiv=1.0
+         if ((itype(j,5).eq.1).or.(itype(j,5).eq.3)) ndiv=2.0
+
+        wh2o=78*ndiv
         wdip =vcatprm(2)
         wdip=wdip/wh2o
         wquad1 =vcatprm(3)