working gradient

diff --git a/source/unres/compare.F90 b/source/unres/compare.F90
index 34f8a2a..8e6c0af 100644 (file)
--- a/source/unres/compare.F90
+++ b/source/unres/compare.F90
@@ -300,6 +300,7 @@
       ncont=0
       ees=0.0
       evdw=0.0
+      print *, "nntt,nct",nnt,nct-2
       do 1 i=nnt,nct-2
         if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) goto 1
         xi=c(1,i)
@@ -492,6 +493,8 @@
       enddo
 
       call elecont(lprint,ncont,icont)
+      print *,"after elecont"
+      if (nres_molec(1).eq.0) return
 
 ! finding parallel beta
 !d      write (iout,*) '------- looking for parallel beta -----------'
@@ -742,7 +745,7 @@
        write(12,'(a37)') "DefPropRes 'coil' '! (helix | sheet)'"
        write(12,'(a20)') "XMacStand ribbon.mac"
          
-        
+      if (nres_molec(1).eq.0) return
        write(iout,*) 'UNRES seq:'
        do j=1,nbfrag
         write(iout,*) 'beta ',(bfrag(i,j),i=1,4)
@@ -815,6 +818,7 @@
 !     &             obr,non_conv)
 !        rms=dsqrt(rms)
         call rmsd(rms)
+        print *,"before contact"
 !elte(iout,*) "rms_nacc before contact"
         call contact(.false.,ncont,icont,co)
         frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
@@ -877,6 +881,7 @@
 !      else
 !      do kkk=1,nperm
       iatom=0
+      print *,nz_start,nz_end,nstart_seq-nstart_sup
       do i=nz_start,nz_end
         iatom=iatom+1
         iti=itype(i,1)

diff --git a/source/unres/control.F90 b/source/unres/control.F90
index d70d144..24ea5e5 100644 (file)
--- a/source/unres/control.F90
+++ b/source/unres/control.F90
@@ -728,7 +728,9 @@
       ispp=4 !?? wham ispp=2
 #ifdef MPI
 ! Now partition the electrostatic-interaction array
-      if (itype(nres_molec(1),1).eq.ntyp1_molec(1)) then
+      if (nres_molec(1).eq.0) then  
+       npept=0
+      elseif (itype(nres_molec(1),1).eq.ntyp1_molec(1)) then
       npept=nres_molec(1)-nnt-1
       else
       npept=nres_molec(1)-nnt
@@ -787,7 +789,7 @@
           ijunk,ielstart_nucl(i),ielend_nucl(i),*113)
       enddo ! i 
   113 continue
-      if (iatel_s.eq.0) iatel_s=1
+      if (iatel_s_nucl.eq.0) iatel_s_nucl=1
 
       nele_int_tot_vdw=(npept-2)*(npept-2+1)/2
 !      write (iout,*) "nele_int_tot_vdw",nele_int_tot_vdw
@@ -810,7 +812,7 @@
       if (iatel_s_vdw.eq.0) iatel_s_vdw=1
    15 continue
       if (iatel_s.eq.0) iatel_s=1
-
+      if (iatel_s_vdw.eq.0) iatel_s_vdw=1
       nele_int_tot_vdw_nucl=(npept_nucl-2)*(npept_nucl-2+1)/2
 !      write (iout,*) "nele_int_tot_vdw",nele_int_tot_vdw
       call int_bounds(nele_int_tot_vdw_nucl,my_ele_inds_vdw_nucl,&
@@ -831,7 +833,7 @@
 !        write (iout,*) i," ielstart_vdw",ielstart_vdw(i),
 !     &   " ielend_vdw",ielend_vdw(i)
       enddo ! i 
-      if (iatel_s_vdw.eq.0) iatel_s_vdw=1
+      if (iatel_s_vdw.eq.0) iatel_s_vdw_nucl=1
   115 continue
 
 #else
@@ -934,6 +936,7 @@
       enddo ! i
   114 continue
       print *, "after inloop3",iatscp_s_nucl,iatscp_e_nucl
+      if (iatscp_s_nucl.eq.0) iatscp_s_nucl=1
 #else
       iatscp_s=nnt
       iatscp_e=nct_molec(1)-1

diff --git a/source/unres/data/energy_data.f90 b/source/unres/data/energy_data.f90
index bbc2ed2..baf7f55 100644 (file)
--- a/source/unres/data/energy_data.f90
+++ b/source/unres/data/energy_data.f90
@@ -66,7 +66,7 @@
        wscloc,wcorr,wcorr4,wcorr5,wcorr6,wsccor,wel_loc,wturn3,wturn4,&
        wturn6,wvdwpp,wliptran,wshield,lipscale,wtube, &
        wbond_nucl,wang_nucl,wcorr_nucl,wcorr3_nucl,welpp,wtor_nucl,&
-       wtor_d_nucl,welsb,wsbloc,wvdwsb,welpsb,wvdwpsb
+       wtor_d_nucl,welsb,wsbloc,wvdwsb,welpsb,wvdwpp_nucl,wvdwpsb
 #ifdef CLUSTER
       real(kind=8) :: scalscp
 #endif

diff --git a/source/unres/energy.f90 b/source/unres/energy.f90
index 881497c..9d92c63 100644 (file)
--- a/source/unres/energy.f90
+++ b/source/unres/energy.f90
@@ -128,7 +128,8 @@
 !-----------------------------NUCLEIC GRADIENT
       real(kind=8),dimension(:,:),allocatable  ::gradb_nucl,gradbx_nucl, &
         gvdwpsb1,gelpp,gvdwpsb,gelsbc,gelsbx,gvdwsbx,gvdwsbc,gsbloc,&
-        gsblocx,gradcorr_nucl,gradxorr_nucl,gradcorr3_nucl,gradxorr3_nucl
+        gsblocx,gradcorr_nucl,gradxorr_nucl,gradcorr3_nucl,gradxorr3_nucl,&
+        gvdwpp_nucl
 !      real(kind=8),dimension(:,:),allocatable :: gloc,gloc_x !(maxvar,2)
       real(kind=8),dimension(:,:),allocatable :: gel_loc,gel_loc_long,&
         gcorr3_turn,gcorr4_turn,gcorr6_turn,gradb,gradbx !(3,maxres)
@@ -274,6 +275,20 @@
           weights_(17)=wbond
           weights_(18)=scal14
           weights_(21)=wsccor
+          weights_(26)=wvdwpp_nucl
+          weights_(27)=welpp
+          weights_(28)=wvdwpsb
+          weights_(29)=welpsb
+          weights_(30)=wvdwsb
+          weights_(31)=welsb
+          weights_(32)=wbond_nucl
+          weights_(33)=wang_nucl
+          weights_(34)=wsbloc
+          weights_(35)=wtor_nucl
+          weights_(36)=wtor_d_nucl
+          weights_(37)=wcorr_nucl
+          weights_(38)=wcorr3_nucl
+
 ! FG Master broadcasts the WEIGHTS_ array
           call MPI_Bcast(weights_(1),n_ene,&
              MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR)
@@ -300,6 +315,20 @@
           wbond=weights(17)
           scal14=weights(18)
           wsccor=weights(21)
+          wvdwpp_nucl =weights(26)
+          welpp  =weights(27)
+          wvdwpsb=weights(28)
+          welpsb =weights(29)
+          wvdwsb =weights(30)
+          welsb  =weights(31)
+          wbond_nucl  =weights(32)
+          wang_nucl   =weights(33)
+          wsbloc =weights(34)
+          wtor_nucl   =weights(35)
+          wtor_d_nucl =weights(36)
+          wcorr_nucl  =weights(37)
+          wcorr3_nucl =weights(38)
+
         endif
         time_Bcast=time_Bcast+MPI_Wtime()-time00
         time_Bcastw=time_Bcastw+MPI_Wtime()-time00
@@ -365,7 +394,7 @@
        if (shield_mode.eq.2) then
                  call set_shield_fac2
        endif
-       print *,"AFTER EGB",ipot,evdw
+!       print *,"AFTER EGB",ipot,evdw
 !mc
 !mc Sep-06: egb takes care of dynamic ss bonds too
 !mc
@@ -431,7 +460,7 @@
 ! Calculate the bond-stretching energy
 !
       call ebond(estr)
-       print *,"EBOND",estr
+!       print *,"EBOND",estr
 !       write(iout,*) "in etotal afer ebond",ipot
 
 ! 
@@ -448,6 +477,7 @@
         call ebend(ebe,ethetacnstr)
       else
         ebe=0
+        ethetacnstr=0
       endif
 !      print *,"Processor",myrank," computed UB"
 !
@@ -549,7 +579,7 @@
        etube=0.0d0
       endif
 !--------------------------------------------------------
-      print *,"before",ees,evdw1,ecorr
+!      print *,"before",ees,evdw1,ecorr
       call ebond_nucl(estr_nucl)
       call ebend_nucl(ebe_nucl)
       call etor_nucl(etors_nucl)
@@ -559,7 +589,7 @@
       call esb(esbloc)
       call multibody_hb_nucl(ecorr_nucl,ecorr3_nucl,n_corr,n_corr1)
 
-      print *,"after ebend", ebe_nucl
+!      print *,"after ebend", ebe_nucl
 #ifdef TIMING
       time_enecalc=time_enecalc+MPI_Wtime()-time00
 #endif
@@ -753,7 +783,7 @@
        +wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran+wtube*etube&
        +Eafmforce+ethetacnstr  &
        +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
-       +wvdwpp*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
+       +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
        +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
        +wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_nucl
 #else
@@ -765,7 +795,7 @@
        +wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran+wtube*etube&
        +Eafmforce+ethetacnstr &
        +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
-       +wvdwpp*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
+       +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
        +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
        +wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_nucl
 #endif
@@ -952,7 +982,7 @@
         edihcnstr,ethetacnstr,ebr*nss,&
         Uconst,eliptran,wliptran,Eafmforce,etube,wtube, & ! till now protein
         estr_nucl,wbond_nucl,ebe_nucl,wang_nucl, &
-        evdwpp,wvdwpp,eespp,welpp,evdwpsb,wvdwpsb,eelpsb,welpsb,&
+        evdwpp,wvdwpp_nucl,eespp,welpp,evdwpsb,wvdwpsb,eelpsb,welpsb,&
         evdwsb,wvdwsb,eelsb,welsb,esbloc,wsbloc,etors_nucl,wtor_nucl,&
         etors_d_nucl,wtor_d_nucl,ecorr_nucl,wcorr_nucl,&
         ecorr3_nucl,wcorr3_nucl, &
@@ -1008,7 +1038,7 @@
         ethetacnstr,ebr*nss,Uconst,eliptran,wliptran,Eafmforc,     &
         etube,wtube, &
         estr_nucl,wbond_nucl, ebe_nucl,wang_nucl,&
-        evdwpp,wvdwpp,eespp,welpp,evdwpsb,wvdwpsb,eelpsb,welpsb&
+        evdwpp,wvdwpp_nucl,eespp,welpp,evdwpsb,wvdwpsb,eelpsb,welpsb&
         evdwsb,wvdwsb,eelsb,welsb,esbloc,wsbloc,etors_nucl,wtor_nucl&
         etors_d_nucl,wtor_d_nucl,ecorr_nucl,wcorr_nucl,&
         ecorr3_nucl,wcorr3_nucl, &
@@ -3067,7 +3097,7 @@
 !
 ! Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
 !
-      print *,"iatel_s,iatel_e,",iatel_s,iatel_e
+!      print *,"iatel_s,iatel_e,",iatel_s,iatel_e
       do i=iatel_s,iatel_e
         if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) cycle
         dxi=dc(1,i)
@@ -5461,6 +5491,8 @@
         if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*E_tc*dthetg2
         gloc(nphi+i-2,icg)=wang*(E_theta+E_tc*dthett)
       enddo
+!      print *,ithetaconstr_start,ithetaconstr_end,"TU"
+
 ! Ufff.... We've done all this!!!
       return
       end subroutine ebend
@@ -5762,7 +5794,7 @@
 !-----------thete constrains
 !      if (tor_mode.ne.2) then
       ethetacnstr=0.0d0
-!C      print *,ithetaconstr_start,ithetaconstr_end,"TU"
+!      print *,ithetaconstr_start,ithetaconstr_end,"TU"
       do i=ithetaconstr_start,ithetaconstr_end
         itheta=itheta_constr(i)
         thetiii=theta(itheta)
@@ -10410,7 +10442,12 @@
                      +wturn4*gshieldc_t4(j,i)&
                      +wel_loc*gshieldc_ll(j,i)&
                      +wtube*gg_tube(j,i) &
-                     +wbond_nucl*gradb_nucl(j,i)
+                     +wvdwpp_nucl*gvdwpp_nucl(j,i)+welpp*gelpp(j,i)+ &
+                     wvdwpsb*(gvdwpsb(j,i)+gvdwpsb1(j,i))+ &
+                     wvdwsb*gvdwsbc(j,i)+welsb*gelsbc(j,i)+ &
+                     wcorr_nucl*gradcorr_nucl(j,i)&
+                     +wcorr3_nucl*gradcorr3_nucl(j,i)
+
         enddo
       enddo 
 #else
@@ -10433,8 +10470,11 @@
                      +wturn4*gshieldc_t4(j,i) &
                      +wel_loc*gshieldc_ll(j,i)&
                      +wtube*gg_tube(j,i) &
-                     +wbond_nucl*gradb_nucl(j,i)
-
+                     +wvdwpp_nucl*gvdwpp_nucl(j,i)+welpp*gelpp(j,i)+ &
+                     wvdwpsb*(gvdwpsb(j,i)+gvdwpsb1(j,i))+ &
+                     wvdwsb*gvdwsbc(j,i)+welsb*gelsbc(j,i)+ &
+                     wcorr_nucl*gradcorr_nucl(j,i)
+                     +wcorr3_nucl*gradcorr3_nucl(j,i)
         enddo
       enddo 
 #endif
@@ -10591,7 +10631,46 @@
                      +wel_loc*gshieldc_ll(j,i) &
                      +wel_loc*gshieldc_loc_ll(j,i) &
                      +wtube*gg_tube(j,i) &
-                     +wbond_nucl*gradb_nucl(j,i)
+                     +0.5d0*(wvdwpp_nucl*gvdwpp_nucl(j,i)+welpp*gelpp(j,i)&
+                     +wvdwpsb*gvdwpsb1(j,i))&
+                     +wbond_nucl*gradb_nucl(j,i)+wsbloc*gsbloc(j,i)
+
+!                 if ((i.le.2).and.(i.ge.1)) print *,gradc(j,i,icg),&
+!                      gradbufc(j,i),welec*gelc(j,i), &
+!                      wel_loc*gel_loc(j,i), &
+!                      wscp*gvdwc_scpp(j,i), &
+!                      welec*gelc_long(j,i),wvdwpp*gvdwpp(j,i), &
+!                      wel_loc*gel_loc_long(j,i), &
+!                      wcorr*gradcorr_long(j,i), &
+!                      wcorr5*gradcorr5_long(j,i), &
+!                      wcorr6*gradcorr6_long(j,i), &
+!                      wturn6*gcorr6_turn_long(j,i), &
+!                      wbond*gradb(j,i), &
+!                      wcorr*gradcorr(j,i), &
+!                      wturn3*gcorr3_turn(j,i), &
+!                      wturn4*gcorr4_turn(j,i), &
+!                      wcorr5*gradcorr5(j,i), &
+!                      wcorr6*gradcorr6(j,i), &
+!                      wturn6*gcorr6_turn(j,i), &
+!                      wsccor*gsccorc(j,i) &
+!                     ,wscloc*gscloc(j,i)  &
+!                     ,wliptran*gliptranc(j,i) &
+!                     ,gradafm(j,i) &
+!                     ,welec*gshieldc(j,i) &
+!                     ,welec*gshieldc_loc(j,i) &
+!                     ,wcorr*gshieldc_ec(j,i) &
+!                     ,wcorr*gshieldc_loc_ec(j,i) &
+!                     ,wturn3*gshieldc_t3(j,i) &
+!                     ,wturn3*gshieldc_loc_t3(j,i) &
+!                     ,wturn4*gshieldc_t4(j,i) &
+!                     ,wturn4*gshieldc_loc_t4(j,i) &
+!                     ,wel_loc*gshieldc_ll(j,i) &
+!                     ,wel_loc*gshieldc_loc_ll(j,i) &
+!                     ,wtube*gg_tube(j,i) &
+!                     ,wbond_nucl*gradb_nucl(j,i) &
+!                     ,wvdwpp_nucl*gvdwpp_nucl(j,i),welpp*gelpp(j,i),&
+!                     wvdwpsb*gvdwpsb1(j,i)&
+!                     ,wbond_nucl*gradb_nucl(j,i),wsbloc*gsbloc(j,i)
 
 
 
@@ -10627,7 +10706,10 @@
                      +wel_loc*gshieldc_ll(j,i) &
                      +wel_loc*gshieldc_loc_ll(j,i) &
                      +wtube*gg_tube(j,i) &
-                     +wbond_nucl*gradb_nucl(j,i) 
+                     +wbond_nucl*gradb_nucl(j,i) &
+                     +0.5d0*(wvdwpp_nucl*gvdwpp_nucl(j,i)+welpp*gelpp(j,i)&
+                     +wvdwpsb*gvdwpsb1(j,i))&
+                     +wsbloc*gsbloc(j,i)
 
 
 
@@ -10645,10 +10727,12 @@
                        +wturn4*gshieldx_t4(j,i) &
                        +wel_loc*gshieldx_ll(j,i)&
                        +wtube*gg_tube_sc(j,i)   &
-                       +wbond_nucl*gradbx_nucl(j,i) 
-
-
-
+                       +wbond_nucl*gradbx_nucl(j,i) &
+                       +wvdwsb*gvdwsbx(j,i) &
+                       +welsb*gelsbx(j,i) &
+                       +wcorr_nucl*gradxorr_nucl(j,i)&
+                       +wcorr3_nucl*gradxorr3_nucl(j,i) &
+                       +wsbloc*gsblocx(j,i)
         enddo
       enddo 
 #ifdef DEBUG
@@ -10675,7 +10759,7 @@
 #ifdef MPI
       if (nfgtasks.gt.1) then
         do j=1,3
-          do i=1,nres
+          do i=0,nres
             gradbufc(j,i)=gradc(j,i,icg)
             gradbufx(j,i)=gradx(j,i,icg)
           enddo
@@ -10704,7 +10788,7 @@
         time00=MPI_Wtime()
         call MPI_Reduce(gradbufc(1,0),gradc(1,0,icg),3*nres+3,&
           MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
-        call MPI_Reduce(gradbufx(1,1),gradx(1,1,icg),3*nres,&
+        call MPI_Reduce(gradbufx(1,0),gradx(1,0,icg),3*nres+3,&
           MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
         call MPI_Reduce(glocbuf(1),gloc(1,icg),4*nres,&
           MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
@@ -10713,6 +10797,7 @@
           MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
         time_reduce=time_reduce+MPI_Wtime()-time00
 !#define DEBUG
+!          print *,"gradbuf",gradbufc(1,1),gradc(1,1,icg)
 #ifdef DEBUG
       write (iout,*) "gloc_sc after reduce"
       do i=1,nres
@@ -11650,7 +11735,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
         enddo
         call zerograd
         call etotal_short(energia)
-!el        call enerprint(energia)
+        call enerprint(energia)
         call flush(iout)
         write (iout,*) "enter cartgrad"
         call flush(iout)
@@ -11857,6 +11942,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
         do i=1,nres
           do j=1,3
             grad_s(j,i)=gcart(j,i)
+!            if (i.le.2) print *,"tu?!",gcart(j,i),grad_s(j,i),gxcart(j,i)
             grad_s(j+3,i)=gxcart(j,i)
           enddo
         enddo
@@ -11934,6 +12020,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
           if (.not.split_ene) then
             call etotal(energia1)
             etot1=energia1(0)
+!            call enerprint(energia1)
           else
 !- split gradient
             call etotal_long(energia1)
@@ -12062,11 +12149,11 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
       call var_to_geom(nvar,x)
       call chainbuild
       icall=1
-      print *,'ICG=',ICG
+!      print *,'ICG=',ICG
       call etotal(energia)
       etot = energia(0)
 !el      call enerprint(energia)
-      print *,'ICG=',ICG
+!      print *,'ICG=',ICG
 #ifdef MPL
       if (MyID.ne.BossID) then
         call mp_bcast(x(1),8*(nvar+3),BossID,fgGroupID)
@@ -16070,6 +16157,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
         do j=1,3
           gcart(j,i)=gradc(j,i,icg)
           gxcart(j,i)=gradx(j,i,icg)
+!          if (i.le.2) print *,"gcart_one",gcart(j,i),gradc(j,i,icg)
         enddo
 #ifdef DEBUG
         write (iout,'(i5,2(3f10.5,5x),f10.5)') i,(gcart(j,i),j=1,3),&
@@ -16080,6 +16168,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
       time01=MPI_Wtime()
 #endif
       call int_to_cart
+!       print *,"gcart_two",gcart(2,2),gradc(2,2,icg)
+
 #ifdef TIMING
       time_inttocart=time_inttocart+MPI_Wtime()-time01
 #endif
@@ -16230,6 +16320,25 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
           gradafm(j,i)=0.0d0
           gradb_nucl(j,i)=0.0d0
           gradbx_nucl(j,i)=0.0d0
+          gvdwpp_nucl(j,i)=0.0d0
+          gvdwpp(j,i)=0.0d0
+          gelpp(j,i)=0.0d0
+          gvdwpsb(j,i)=0.0d0
+          gvdwpsb1(j,i)=0.0d0
+          gvdwsbc(j,i)=0.0d0
+          gvdwsbx(j,i)=0.0d0
+          gelsbc(j,i)=0.0d0
+          gradcorr_nucl(j,i)=0.0d0
+          gradcorr3_nucl(j,i)=0.0d0
+          gradxorr_nucl(j,i)=0.0d0
+          gradxorr3_nucl(j,i)=0.0d0
+          gelsbx(j,i)=0.0d0
+          gsbloc(j,i)=0.0d0
+          gsblocx(j,i)=0.0d0
+        enddo
+       enddo
+      do i=0,nres
+        do j=1,3
           do intertyp=1,3
            gloc_sc(intertyp,i,icg)=0.0d0
           enddo
@@ -19604,6 +19713,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
       allocate(gradxorr_nucl(3,-1:nres))
       allocate(gradcorr3_nucl(3,-1:nres))
       allocate(gradxorr3_nucl(3,-1:nres))
+      allocate(gvdwpp_nucl(3,-1:nres))
 
 !(3,maxres)
       allocate(grad_shield_side(3,50,nres))
@@ -19813,14 +19923,14 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
           if(energy_dec)write(iout,*) "estr_nucl_bb" , i,vbld(i),&
           vbldp0_nucl,diff,AKP_nucl*diff*diff
           estr_nucl=estr_nucl+diff*diff
-          print *,estr_nucl
+!          print *,estr_nucl
           do j=1,3
             gradb_nucl(j,i-1)=AKP_nucl*diff*dc(j,i-1)/vbld(i)
           enddo
 !c          write (iout,'(i5,3f10.5)') i,(gradb(j,i-1),j=1,3)
       enddo
       estr_nucl=0.5d0*AKP_nucl*estr_nucl
-      print *,"partial sum", estr_nucl,AKP_nucl
+!      print *,"partial sum", estr_nucl,AKP_nucl
 
       if (energy_dec) &
       write (iout,*) "ibondp_start,ibondp_end",&
@@ -19839,7 +19949,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
            write (iout,*) "estr_nucl_sc", i,iti,vbld(i+nres),vbldsc0_nucl(1,iti),diff, &
            AKSC_nucl(1,iti),AKSC_nucl(1,iti)*diff*diff
             estr_nucl=estr_nucl+0.5d0*AKSC_nucl(1,iti)*diff*diff
-            print *,estr_nucl
+!            print *,estr_nucl
             do j=1,3
               gradbx_nucl(j,i)=AKSC_nucl(1,iti)*diff*dc(j,i+nres)/vbld(i+nres)
             enddo
@@ -19869,7 +19979,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
             enddo
             estr_nucl=estr_nucl+uprod/usum
             do j=1,3
-             gradbx(j,i)=usumsqder/(usum*usum)*dc(j,i+nres)/vbld(i+nres)
+             gradbx_nucl(j,i)=usumsqder/(usum*usum)*dc(j,i+nres)/vbld(i+nres)
             enddo
         endif
       enddo
@@ -19882,7 +19992,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
       real(kind=8),dimension(nntheterm_nucl+1) :: coskt,sinkt !mmaxtheterm
       real(kind=8),dimension(nsingle_nucl+1) :: cosph1,sinph1,cosph2,sinph2 !maxsingle
       real(kind=8),dimension(ndouble_nucl+1,ndouble_nucl+1) :: cosph1ph2,sinph1ph2 !maxdouble,maxdouble
-      logical :: lprn=.true., lprn1=.false.
+      logical :: lprn=.false., lprn1=.false.
 !el local variables
       integer :: i,k,iblock,ityp1,ityp2,ityp3,l,m
       real(kind=8) :: dethetai,dephii,dephii1,theti2,phii,phii1,ethetai
@@ -20166,7 +20276,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
 !c
 !c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
 !c
-      print *,"iatel_s_nucl,iatel_e_nucl",iatel_s_nucl,iatel_e_nucl
+!      print *,"iatel_s_nucl,iatel_e_nucl",iatel_s_nucl,iatel_e_nucl
       do i=iatel_s_nucl,iatel_e_nucl
         if (itype(i,2).eq.ntyp1_molec(2) .or. itype(i+1,2).eq.ntyp1_molec(2)) cycle
         dxi=dc(1,i)
@@ -20253,8 +20363,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
           ggg(2)=fac*yj
           ggg(3)=fac*zj
           do k=1,3
-            gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
-            gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
+            gvdwpp_nucl(k,i)=gvdwpp_nucl(k,i)-ggg(k)
+            gvdwpp_nucl(k,j)=gvdwpp_nucl(k,j)+ggg(k)
           enddo
 !c phoshate-phosphate electrostatic interactions
           rij=dsqrt(rij)
@@ -20282,7 +20392,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
       do i=nnt,nct
 !c        write (2,*) "i",i," gelpp",(gelpp(k,i),k=1,3)
         do k=1,3
-          gvdwpp(k,i)=6*gvdwpp(k,i)
+          gvdwpp_nucl(k,i)=6*gvdwpp_nucl(k,i)
 !c          gelpp(k,i)=332.0d0*gelpp(k,i)
           gelpp(k,i)=AEES*gelpp(k,i)
         enddo
@@ -20311,7 +20421,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
 !cd    write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
       eelpsb=0.0d0
       evdwpsb=0.0d0
-      print *,"iatscp_s_nucl,iatscp_e_nucl",iatscp_s_nucl,iatscp_e_nucl
+!      print *,"iatscp_s_nucl,iatscp_e_nucl",iatscp_s_nucl,iatscp_e_nucl
       do i=iatscp_s_nucl,iatscp_e_nucl
         if (itype(i,2).eq.ntyp1_molec(2) &
          .or. itype(i+1,2).eq.ntyp1_molec(2)) cycle
@@ -20910,8 +21020,9 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
         enddo
         do j = 1,3
           z_prime(j) = -uz(j,i-1)
+!           z_prime(j)=0.0
         enddo
-
+       
         xx=0.0d0
         yy=0.0d0
         zz=0.0d0
@@ -20941,8 +21052,10 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
 #endif
         sumene = enesc_nucl(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1))
         esbloc = esbloc + sumene
-        if (energy_dec) write(iout,*) "i",i," esbloc",sumene,esbloc,xx,yy,zz
-        if (energy_dec) write(iout,*) "x",(x(k),k=1,9)
+        sumene2= enesc_nucl(x,xx,yy,0.0d0,cost2tab(i+1),sint2tab(i+1))
+!        print *,"enecomp",sumene,sumene2
+!        if (energy_dec) write(iout,*) "i",i," esbloc",sumene,esbloc,xx,yy,zz
+!        if (energy_dec) write(iout,*) "x",(x(k),k=1,9)
 #ifdef DEBUG
         write (2,*) "x",(x(k),k=1,9)
 !C
@@ -21054,12 +21167,13 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
 !c     &    dt_dci(k)
 !c         write (iout,*) "k",k," dxx_XYZ",dxx_XYZ(k)," dyy_XYZ",
 !c     &    dyy_XYZ(k)," dzz_XYZ",dzz_XYZ(k) 
-         gsbloc(k,i-1)=gsbloc(k,i-1)+de_dxx*dxx_ci1(k) &
-         +de_dyy*dyy_ci1(k)+de_dzz*dzz_ci1(k)+de_dt*dt_dCi1(k)
-         gsbloc(k,i)=gsbloc(k,i)+de_dxx*dxx_Ci(k) &
-         +de_dyy*dyy_Ci(k)+de_dzz*dzz_Ci(k)+de_dt*dt_dCi(k)
+         gsbloc(k,i-1)=gsbloc(k,i-1)+(de_dxx*dxx_ci1(k) &
+         +de_dyy*dyy_ci1(k)+de_dzz*dzz_ci1(k)+de_dt*dt_dCi1(k))
+         gsbloc(k,i)=gsbloc(k,i)+(de_dxx*dxx_Ci(k) &
+         +de_dyy*dyy_Ci(k)+de_dzz*dzz_Ci(k)+de_dt*dt_dCi(k))
          gsblocx(k,i)=                 de_dxx*dxx_XYZ(k)&
          +de_dyy*dyy_XYZ(k)+de_dzz*dzz_XYZ(k)
+!         print *,i,de_dxx*dxx_ci1(k)+de_dyy*dyy_ci1(k),de_dzz*dzz_ci1(k)*2
        enddo
 !c       write(iout,*) "ENERGY GRAD = ", (gsbloc(k,i-1),k=1,3),
 !c     &  (gsbloc(k,i),k=1,3),(gsblocx(k,i),k=1,3)  
@@ -21233,21 +21347,21 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
       ecorr=0.0D0
       ecorr3=0.0d0
 !C Remove the loop below after debugging !!!
-      do i=nnt_molec(2),nct_molec(2)
-        do j=1,3
-          gradcorr_nucl(j,i)=0.0D0
-          gradxorr_nucl(j,i)=0.0D0
-          gradcorr3_nucl(j,i)=0.0D0
-          gradxorr3_nucl(j,i)=0.0D0
-        enddo
-      enddo
-      print *,"iatsc_s_nucl,iatsc_e_nucl",iatsc_s_nucl,iatsc_e_nucl
+!      do i=nnt_molec(2),nct_molec(2)
+!        do j=1,3
+!          gradcorr_nucl(j,i)=0.0D0
+!          gradxorr_nucl(j,i)=0.0D0
+!          gradcorr3_nucl(j,i)=0.0D0
+!          gradxorr3_nucl(j,i)=0.0D0
+!        enddo
+!      enddo
+!      print *,"iatsc_s_nucl,iatsc_e_nucl",iatsc_s_nucl,iatsc_e_nucl
 !C Calculate the local-electrostatic correlation terms
       do i=iatsc_s_nucl,iatsc_e_nucl
         i1=i+1
         num_conti=num_cont_hb(i)
         num_conti1=num_cont_hb(i+1)
-        print *,i,num_conti,num_conti1
+!        print *,i,num_conti,num_conti1
         do jj=1,num_conti
           j=jcont_hb(jj,i)
           do kk=1,num_conti1

diff --git a/source/unres/geometry.f90 b/source/unres/geometry.f90
index e0b37d3..653e595 100644 (file)
--- a/source/unres/geometry.f90
+++ b/source/unres/geometry.f90
@@ -3240,6 +3240,8 @@
 !   calculating dE/ddc1  
 !el local variables
        integer :: j,i
+!       print *,"gloc",gloc(:,:)
+!       print *, "gcart",gcart(:,:)
        if (nres.lt.3) go to 18
        do j=1,3
          gcart(j,1)=gcart(j,1)+gloc(1,icg)*dphi(j,1,4) &
@@ -3331,7 +3333,8 @@
         enddo
 !   The side-chain vector derivatives
         do i=2,nres-1
-         if(itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then   
+         if(itype(i,1).ne.10 .and.  &
+           itype(i,molnum(i)).ne.ntyp1_molec(molnum(i))) then  
             do j=1,3   
               gxcart(j,i)=gxcart(j,i)+gloc(ialph(i,1),icg)*dalpha(j,3,i) &
               +gloc(ialph(i,1)+nside,icg)*domega(j,3,i)
@@ -3350,7 +3353,8 @@
 !       enddo
        if (nres.lt.2) return
        if ((nres.lt.3).and.(itype(1,1).eq.10)) return
-       if ((itype(1,1).ne.10).and.(itype(1,1).ne.ntyp1)) then
+       if ((itype(1,1).ne.10).and. &
+        (itype(1,molnum(1)).ne.ntyp1_molec(molnum(1)))) then
         do j=1,3
 !c Derviative was calculated for oposite vector of side chain therefore
 ! there is "-" sign before gloc_sc
@@ -3358,7 +3362,8 @@
            dtauangle(j,1,1,3)
          gcart(j,1)=gcart(j,1)+gloc_sc(1,0,icg)* &
            dtauangle(j,1,2,3)
-          if ((itype(2,1).ne.10).and.(itype(2,1).ne.ntyp1)) then
+          if ((itype(2,1).ne.10).and. &
+        (itype(2,molnum(2)).ne.ntyp1_molec(molnum(2)))) then
          gxcart(j,1)= gxcart(j,1) &
                      -gloc_sc(3,0,icg)*dtauangle(j,3,1,3)
          gcart(j,1)=gcart(j,1)+gloc_sc(3,0,icg)* &
@@ -3366,7 +3371,8 @@
           endif
        enddo
        endif
-         if ((nres.ge.3).and.(itype(3,1).ne.10).and.(itype(3,1).ne.ntyp1)) &
+         if ((nres.ge.3).and.(itype(3,molnum(3)).ne.10).and.&
+         (itype(3,molnum(3)).ne.ntyp1_molec(molnum(3)))) &
       then
          do j=1,3
          gcart(j,1)=gcart(j,1)+gloc_sc(2,1,icg)*dtauangle(j,2,1,4)
@@ -3378,10 +3384,11 @@
 
 !     Calculating the remainder of dE/ddc2
        do j=1,3
-         if((itype(2,1).ne.10).and.(itype(2,1).ne.ntyp1)) then
+         if((itype(2,1).ne.10).and. &
+           (itype(2,molnum(2)).ne.ntyp1_molec(molnum(2)))) then
            if (itype(1,1).ne.10) gxcart(j,2)=gxcart(j,2)+ &
                                gloc_sc(3,0,icg)*dtauangle(j,3,3,3)
-        if ((itype(3,1).ne.10).and.(nres.ge.3).and.(itype(3,1).ne.ntyp1)) &
+        if ((itype(3,1).ne.10).and.(nres.ge.3).and.(itype(3,molnum(3)).ne.ntyp1_molec(3))) &
          then
            gxcart(j,2)=gxcart(j,2)-gloc_sc(3,1,icg)*dtauangle(j,3,1,4)
 !c                  the   - above is due to different vector direction
@@ -3395,7 +3402,8 @@
 !           write(iout,*) gloc_sc(1,1,icg),dtauangle(j,1,1,4),"gx"
           endif
          endif
-         if ((itype(1,1).ne.10).and.(itype(1,1).ne.ntyp1)) then
+         if ((itype(1,1).ne.10).and.&
+         (itype(1,molnum(1)).ne.ntyp1_molec(molnum(1)))) then
           gcart(j,2)=gcart(j,2)+gloc_sc(1,0,icg)*dtauangle(j,1,3,3)
 !           write(iout,*)  gloc_sc(1,0,icg),dtauangle(j,1,3,3)
         endif
@@ -3415,7 +3423,8 @@
         do i=3,nres-2
          do j=1,3
 !          write(iout,*) "before", gcart(j,i)
-          if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
+          if ((itype(i,1).ne.10).and.&
+          (itype(i,molnum(i)).ne.ntyp1_molec(molnum(i)))) then
           gxcart(j,i)=gxcart(j,i)+gloc_sc(2,i-2,icg) &
           *dtauangle(j,2,3,i+1) &
           -gloc_sc(1,i-1,icg)*dtauangle(j,1,1,i+2)
@@ -3454,7 +3463,8 @@
 !  Setting dE/ddnres-1       
       if(nres.ge.4) then
          do j=1,3
-         if ((itype(nres-1,1).ne.10).and.(itype(nres-1,1).ne.ntyp1)) then
+         if ((itype(nres-1,1).ne.10).and.&
+       (itype(nres-1,molnum(nres-1)).ne.ntyp1_molec(molnum(nres-1)))) then
          gxcart(j,nres-1)=gxcart(j,nres-1)+gloc_sc(2,nres-3,icg) &
           *dtauangle(j,2,3,nres)
 !          write (iout,*) "gxcart(nres-1)", gloc_sc(2,nres-3,icg),
@@ -3463,18 +3473,20 @@
         gxcart(j,nres-1)=gxcart(j,nres-1)+gloc_sc(3,nres-3,icg) &
           *dtauangle(j,3,3,nres)
           endif
-         if ((itype(nres,1).ne.10).and.(itype(nres,1).ne.ntyp1)) then
+         if ((itype(nres,1).ne.10).and.&
+         (itype(nres,molnum(nres)).ne.ntyp1_molec(molnum(nres)))) then
         gxcart(j,nres-1)=gxcart(j,nres-1)-gloc_sc(3,nres-2,icg) &
           *dtauangle(j,3,1,nres+1)
         gcart(j,nres-1)=gcart(j,nres-1)+gloc_sc(3,nres-2,icg) &
           *dtauangle(j,3,2,nres+1)
           endif
          endif
-         if ((itype(nres-2,1).ne.10).and.(itype(nres-2,1).ne.ntyp1)) then
+         if ((itype(nres-2,1).ne.10).and.&
+         (itype(nres-2,molnum(nres-2)).ne.ntyp1_molec(molnum(nres-2)))) then
             gcart(j,nres-1)=gcart(j,nres-1)+gloc_sc(1,nres-3,icg)* &
          dtauangle(j,1,3,nres)
          endif
-          if ((itype(nres,1).ne.10).and.(itype(nres,1).ne.ntyp1)) then
+          if ((itype(nres,1).ne.10).and.(itype(nres,molnum(nres)).ne.ntyp1_molec(molnum(nres)))) then
             gcart(j,nres-1)=gcart(j,nres-1)+gloc_sc(2,nres-2,icg)* &
            dtauangle(j,2,2,nres+1)
 !           write (iout,*) "gcart(nres-1)", gloc_sc(2,nres-2,icg),
@@ -3484,7 +3496,7 @@
       endif
 !  Settind dE/ddnres       
        if ((nres.ge.3).and.(itype(nres,1).ne.10).and. &
-          (itype(nres,1).ne.ntyp1))then
+          (itype(nres,molnum(nres)).ne.ntyp1_molec(molnum(nres))))then
        do j=1,3
         gxcart(j,nres)=gxcart(j,nres)+gloc_sc(3,nres-2,icg) &
        *dtauangle(j,3,3,nres+1)+gloc_sc(2,nres-2,icg) &
@@ -3492,6 +3504,7 @@
         enddo
        endif
 !   The side-chain vector derivatives
+!       print *,"gcart",gcart(:,:)
       return
       end subroutine int_to_cart
 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)

diff --git a/source/unres/io.f90 b/source/unres/io.f90
index ee588d0..a86fc0b 100644 (file)
--- a/source/unres/io.f90
+++ b/source/unres/io.f90
@@ -720,7 +720,7 @@
 !      integer :: rescode
 !      double precision x(maxvar)
       character(len=256) :: pdbfile
-      character(len=480) :: weightcard
+      character(len=800) :: weightcard
       character(len=80) :: weightcard_t!,ucase
 !      integer,dimension(:),allocatable :: itype_pdb   !(maxres)
 !      common /pizda/ itype_pdb
@@ -769,6 +769,19 @@
       call reada(weightcard,'WTURN6',wturn6,1.0D0)
       call reada(weightcard,'WSCCOR',wsccor,1.0D0)
       call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
+      call reada(weightcard,'WVDWPP_NUCL',wvdwpp_nucl,1.0D0)
+      call reada(weightcard,'WELPP',welpp,1.0d0)
+      call reada(weightcard,'WVDWPSB',wvdwpsb,1.0d0)
+      call reada(weightcard,'WELPSB',welpsb,1.0D0)
+      call reada(weightcard,'WVDWSB',wvdwsb,1.0d0)
+      call reada(weightcard,'WELSB',welsb,1.0D0)
+      call reada(weightcard,'WBOND_NUCL',wbond_nucl,1.0D0)
+      call reada(weightcard,'WANG_NUCL',wang_nucl,1.0D0)
+      call reada(weightcard,'WSBLOC',wsbloc,1.0D0)
+      call reada(weightcard,'WTOR_NUCL',wtor_nucl,1.0D0)
+      call reada(weightcard,'WTORD_NUCL',wtor_d_nucl,1.0D0)
+      call reada(weightcard,'WCORR_NUCL',wcorr_nucl,1.0D0)
+      call reada(weightcard,'WCORR3_NUCL',wcorr3_nucl,1.0D0)
       call reada(weightcard,'WBOND',wbond,1.0D0)
       call reada(weightcard,'WTOR',wtor,1.0D0)
       call reada(weightcard,'WTORD',wtor_d,1.0D0)
@@ -1188,8 +1201,8 @@
 #endif
       nct=nres
       print *,'NNT=',NNT,' NCT=',NCT
-      if (itype(1,1).eq.ntyp1) nnt=2
-      if (itype(nres,1).eq.ntyp1) nct=nct-1
+      if (itype(1,molnum(1)).eq.ntyp1_molec(molnum(1))) nnt=2
+      if (itype(nres,molnum(nres)).eq.ntyp1_molec(molnum(nres))) nct=nct-1
       if (pdbref) then
         if(me.eq.king.or..not.out1file) &
          write (iout,'(a,i3)') 'nsup=',nsup

diff --git a/source/unres/unres.f90 b/source/unres/unres.f90
index 40c0140..e84e1d0 100644 (file)
--- a/source/unres/unres.f90
+++ b/source/unres/unres.f90
@@ -346,7 +346,7 @@
           write (iout,'(a,i20)') '# of energy evaluations:',nfun+1
           write (iout,'(a,f16.3)')'# of energy evaluations/sec:',evals
       else
-        print *,'refstr=',refstr
+        print *,'refstr=',refstr,frac,frac_nn,co
         if (refstr) call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
         print *,"after rms_nac_ncc"
         call briefout(0,etot)