correcation in ion param and dyn_ss
[unres4.git] / source / unres / energy.F90
index 63a76fb..b8587fc 100644 (file)
         if (nfgtasks.gt.1) then 
         call MPI_Bcast(itime_mat,1,MPI_INT,king,FG_COMM,IERROR)
         endif
+       if (nres_molec(1).gt.0) then
        if (mod(itime_mat,imatupdate).eq.0) call make_SCp_inter_list
 !       write (iout,*) "after make_SCp_inter_list"
        if (mod(itime_mat,imatupdate).eq.0) call make_SCSC_inter_list
 !       write (iout,*) "after make_SCSC_inter_list"
 
        if (mod(itime_mat,imatupdate).eq.0) call make_pp_inter_list
+       endif
 !       write (iout,*) "after make_pp_inter_list"
 
 !      print *,'Processor',myrank,' calling etotal ipot=',ipot
       dGCLdOM1=0.0d0 
       dPOLdOM1=0.0d0
 !             write (iout,*) "RWA", g_listscsc_start,g_listscsc_end,i,j
-
+      if (nres_molec(1).eq.0) return
       do icont=g_listscsc_start,g_listscsc_end
       i=newcontlisti(icont)
       j=newcontlistj(icont)
                               'evdw',i,j,evdwij,' ss'
 !              if (energy_dec) write (iout,*) &
 !                              'evdw',i,j,evdwij,' ss'
-             do k=j+1,iend(i,iint)
+             do k=j+1,nres
 !C search over all next residues
               if (dyn_ss_mask(k)) then
 !C check if they are cysteins
       eel_loc=0.0d0 
       eello_turn3=0.0d0
       eello_turn4=0.0d0
+      if (nres_molec(1).eq.0) return
 !
 
       if (icheckgrad.eq.1) then
 #ifdef NEWCORR
         gloc(nphi+i,icg)=gloc(nphi+i,icg)&
                        -(gs13+gsE13+gsEE1)*wturn4&
-       *fac_shield(i)*fac_shield(j)
+       *fac_shield(i)*fac_shield(j) &
+       *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
         gloc(nphi+i+1,icg)= gloc(nphi+i+1,icg)&
                          -(gs23+gs21+gsEE2)*wturn4&
-       *fac_shield(i)*fac_shield(j)
+       *fac_shield(i)*fac_shield(j)&
+       *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
 
         gloc(nphi+i+2,icg)= gloc(nphi+i+2,icg)&
                          -(gs32+gsE31+gsEE3)*wturn4&
-       *fac_shield(i)*fac_shield(j)
+       *fac_shield(i)*fac_shield(j)&
+       *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
 
 !c         gloc(nphi+i+1,icg)=gloc(nphi+i+1,icg)-
 !c     &   gs2
 !d    print '(a)','Enter ESCP'
 !d    write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
 !      do i=iatscp_s,iatscp_e
+      if (nres_molec(1).eq.0) return
        do icont=g_listscp_start,g_listscp_end
         i=newcontlistscpi(icont)
         j=newcontlistscpj(icont)
@@ -19040,14 +19049,12 @@ chip1=chip(itypi)
         if (idssb(i).eq.newihpb(j) .and. &
              jdssb(i).eq.newjhpb(j)) found=.true.
       enddo
-#ifndef CLUST
-#ifndef WHAM
+#if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
 !        write(iout,*) "found",found,i,j
       if (.not.found.and.fg_rank.eq.0) &
           write(iout,'(a15,f12.2,f8.1,2i5)') &
            "SSBOND_BREAK",totT,t_bath,idssb(i),jdssb(i)
 #endif
-#endif
       enddo
 
       do i=1,newnss
@@ -19057,21 +19064,22 @@ chip1=chip(itypi)
         if (newihpb(i).eq.idssb(j) .and. &
              newjhpb(i).eq.jdssb(j)) found=.true.
       enddo
-#ifndef CLUST
-#ifndef WHAM
+#if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
 !        write(iout,*) "found",found,i,j
       if (.not.found.and.fg_rank.eq.0) &
           write(iout,'(a15,f12.2,f8.1,2i5)') &
            "SSBOND_FORM",totT,t_bath,newihpb(i),newjhpb(i)
 #endif
-#endif
       enddo
-
+!#if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
       nss=newnss
       do i=1,nss
       idssb(i)=newihpb(i)
       jdssb(i)=newjhpb(i)
       enddo
+!#else
+!      nss=0
+!#endif
 
       return
       end subroutine dyn_set_nss
@@ -20551,10 +20559,10 @@ chip1=chip(itypi)
 !(3,3,2,maxres)
 ! allocateion of lists JPRDLA
       allocate(newcontlistppi(300*nres))
-      allocate(newcontlistscpi(300*nres))
+      allocate(newcontlistscpi(350*nres))
       allocate(newcontlisti(300*nres))
       allocate(newcontlistppj(300*nres))
-      allocate(newcontlistscpj(300*nres))
+      allocate(newcontlistscpj(350*nres))
       allocate(newcontlistj(300*nres))
 
       return
@@ -22402,7 +22410,7 @@ chip1=chip(itypi)
       enddo
       call to_box(chead(1,1),chead(2,1),chead(3,1))
       call to_box(chead(1,2),chead(2,2),chead(3,2))
-
+!      write(iout,*) "TEST",chead(1,1),chead(2,1),chead(3,1),dc_norm(k, i+nres),d1 
 ! distance 
 !        Rsc_distance(k) = dabs(c(k, i+nres) - c(k, j+nres))
 !         Rsc(k) = Rsc_distance(k) * Rsc_distance(k)
@@ -22455,6 +22463,16 @@ chip1=chip(itypi)
         rij_shift = Rtail - sig + sig0ij
         IF (rij_shift.le.0.0D0) THEN
          evdw = 1.0D20
+      if (evdw.gt.1.0d6) then
+      write (*,'(2(1x,a3,i3),7f7.2)') &
+      restyp(itype(i,1),1),i,restyp(itype(j,1),1),j,&
+      1.0d0/rij,Rtail,Rhead,rij_shift, sig, sig0ij,sigsq
+      write(*,*) facsig,faceps1_inv,om1,chiom1,chi1
+     write(*,*) "ANISO?!",chi1
+!evdwij,Fcav,Ecl,Egb,Epol,Fisocav,Elj,&
+!      Equad,evdwij+Fcav+eheadtail,evdw
+      endif
+
          RETURN
         END IF
         sigder = -sig * sigsq
@@ -22610,6 +22628,12 @@ chip1=chip(itypi)
 !           CALL energy_quad(istate,eheadtail,Ecl,Egb,Epol,Fisocav,Elj,Equad)
        END IF  ! this endif ends the "catch the gly-gly" at the beggining of Fcav
       evdw = evdw  + Fcav + eheadtail
+!      if (evdw.gt.1.0d6) then
+!      write (*,'(2(1x,a3,i3),3f6.2,10f16.7)') &
+!      restyp(itype(i,1),1),i,restyp(itype(j,1),1),j,&
+!      1.0d0/rij,Rtail,Rhead,evdwij,Fcav,Ecl,Egb,Epol,Fisocav,Elj,&
+!      Equad,evdwij+Fcav+eheadtail,evdw
+!      endif
 
        IF (energy_dec) write (iout,'(2(1x,a3,i3),3f6.2,10f16.7)') &
       restyp(itype(i,1),1),i,restyp(itype(j,1),1),j,&
@@ -22790,6 +22814,13 @@ chip1=chip(itypi)
         rij_shift = Rtail - sig + sig0ij
         IF (rij_shift.le.0.0D0) THEN
          evdw = 1.0D20
+!      if (evdw.gt.1.0d6) then
+!      write (*,'(2(1x,a3,i3),6f6.2)') &
+!      restyp(itype(i,1),1),i,restyp(itype(j,1),1),j,&
+!      1.0d0/rij,Rtail,Rhead,rij_shift, sig, sig0ij
+!evdwij,Fcav,Ecl,Egb,Epol,Fisocav,Elj,&
+!      Equad,evdwij+Fcav+eheadtail,evdw
+!      endif
          RETURN
         END IF
         sigder = -sig * sigsq
@@ -22898,7 +22929,12 @@ chip1=chip(itypi)
 !           eheadtail = 0.0d0
 
       evdw = evdw  + Fcav + eheadtail
-
+!      if (evdw.gt.1.0d6) then
+!      write (*,'(2(1x,a3,i3),3f6.2,10f16.7)') &
+!      restyp(itype(i,1),1),i,restyp(itype(j,1),1),j,&
+!      1.0d0/rij,Rtail,Rhead,evdwij,Fcav,Ecl,Egb,Epol,Fisocav,Elj,&
+!      Equad,evdwij+Fcav+eheadtail,evdw
+!      endif
        IF (energy_dec) write (iout,'(2(1x,a3,i3),3f6.2,10f16.7)') &
       restyp(itype(i,1),1),i,restyp(itype(j,1),1),j,&
       1.0d0/rij,Rtail,Rhead,evdwij,Fcav,Ecl,Egb,Epol,Fisocav,Elj,&
@@ -27584,8 +27620,10 @@ chip1=chip(itypi)
            zi=c(3,nres+i)
           call to_box(xi,yi,zi)
            do iint=1,nint_gr(i)
+!           print *,"is it wrong", iint,i
             do j=istart(i,iint),iend(i,iint)
              itypj=iabs(itype(j,1))
+             if (energy_dec) write(iout,*) "LISTA ZAKRES",istart(i,iint),iend(i,iint),iatsc_s,iatsc_e
              if (itypj.eq.ntyp1) cycle
              xj=c(1,nres+j)
              yj=c(2,nres+j)
@@ -27672,7 +27710,7 @@ chip1=chip(itypi)
       include 'mpif.h'
       real*8 :: xi,yi,zi,xj,yj,zj,xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp
       real*8 :: dist_init, dist_temp,r_buff_list
-      integer:: contlistscpi(250*nres),contlistscpj(250*nres)
+      integer:: contlistscpi(350*nres),contlistscpj(350*nres)
 !      integer :: newcontlistscpi(200*nres),newcontlistscpj(200*nres)
       integer i,j,itypi,itypj,subchap,xshift,yshift,zshift,iint,ilist_scp,g_ilist_scp
       integer displ(0:nprocs),i_ilist_scp(0:nprocs),ierr