corrections in wham and clust
authorAdam Sieradzan <adasko@piasek4.chem.univ.gda.pl>
Wed, 21 Oct 2015 09:41:35 +0000 (11:41 +0200)
committerAdam Sieradzan <adasko@piasek4.chem.univ.gda.pl>
Wed, 21 Oct 2015 09:41:35 +0000 (11:41 +0200)
source/cluster/wham/src-M/energy_p_new.F
source/cluster/wham/src-M/main_clust.F
source/cluster/wham/src-M/readrtns.F
source/unres/src_MD-M/energy_p_new_barrier.F
source/wham/src-M/enecalc1.F
source/wham/src-M/energy_p_new.F

index f640679..f2325c6 100644 (file)
@@ -1907,15 +1907,15 @@ cd      write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e
         gcorr_loc(i)=0.0d0
       enddo
       do i=iatel_s,iatel_e
-          if (i.eq.1) then
-           if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1
-     &  .or. itype(i+2).eq.ntyp1) cycle
-          else
+          if (i.eq.1) cycle
+C           if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1
+C     &  .or. itype(i+2).eq.ntyp1) cycle
+C          else
         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
      &  .or. itype(i+2).eq.ntyp1
      &  .or. itype(i-1).eq.ntyp1
      &) cycle
-         endif
+C         endif
         if (itel(i).eq.0) goto 1215
         dxi=dc(1,i)
         dyi=dc(2,i)
@@ -1935,16 +1935,16 @@ cd      write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e
         num_conti=0
 c        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
         do j=ielstart(i),ielend(i)
-          if (j.eq.1) then
-           if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1
-     & .or.itype(j+2).eq.ntyp1
-     &) cycle
-          else
+          if (j.eq.1) cycle
+C           if (itype(j).eq.ntyp1 .or. itype(j+1).eq.Cntyp1
+C     & .or.itype(j+2).eq.ntyp1
+C     &) cycle
+C          else
           if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1
      & .or.itype(j+2).eq.ntyp1
      & .or.itype(j-1).eq.ntyp1
      &) cycle
-         endif
+C         endif
           if (itel(j).eq.0) goto 1216
           ind=ind+1
           iteli=itel(i)
@@ -1975,7 +1975,7 @@ C End diagnostics
           if (yj.lt.0) yj=yj+boxysize
           zj=mod(zj,boxzsize)
           if (zj.lt.0) zj=zj+boxzsize
-      dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+      dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
       xj_safe=xj
       yj_safe=yj
       zj_safe=zj
@@ -1986,7 +1986,7 @@ C End diagnostics
           xj=xj_safe+xshift*boxxsize
           yj=yj_safe+yshift*boxysize
           zj=zj_safe+zshift*boxzsize
-          dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+          dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
           if(dist_temp.lt.dist_init) then
             dist_init=dist_temp
             xj_temp=xj
index 892a6c7..b230c01 100644 (file)
@@ -63,6 +63,7 @@ C
       call openunits
       call parmread
       call read_control
+      write(iout,*) "afert readcontrol"
       call molread
 c      if (refstr) call read_ref_structure(*30)
       do i=1,nres
index ea51dbf..6769fb7 100644 (file)
@@ -118,8 +118,10 @@ C Body
 C
 C Read weights of the subsequent energy terms.
       call card_concat(weightcard)
+      write(iout,*) weightcard
       call reada(weightcard,'WSC',wsc,1.0d0)
       call reada(weightcard,'WLONG',wsc,wsc)
+      write(iout,*) "WLONG=",wsc
       call reada(weightcard,'WSCP',wscp,1.0d0)
       call reada(weightcard,'WELEC',welec,1.0D0)
       call reada(weightcard,'WVDWPP',wvdwpp,welec)
@@ -150,10 +152,12 @@ C Read weights of the subsequent energy terms.
       call reada(weightcard,"V2SS",v2ss,7.61d0)
       call reada(weightcard,"V3SS",v3ss,13.7d0)
       call reada(weightcard,"EBR",ebr,-5.50D0)
+      write (iout,*) "atriss",atriss
       call reada(weightcard,"ATRISS",atriss,0.301D0)
       call reada(weightcard,"BTRISS",btriss,0.021D0)
       call reada(weightcard,"CTRISS",ctriss,1.001D0)
       call reada(weightcard,"DTRISS",dtriss,1.001D0)
+      write(iout,*) "after"
       write (iout,*) "ATRISS=", atriss
       write (iout,*) "BTRISS=", btriss
       write (iout,*) "CTRISS=", ctriss
@@ -213,7 +217,7 @@ C 12/1/95 Added weight for the multi-body term WCORR
      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,
      &  wturn4,wturn6,wsccor
    10 format (/'Energy-term weights (unscaled):'//
-     & 'WSCC=   ',f10.6,' (SC-SC)'/
+     & 'WSC=    ',f10.6,' (SC-SC)'/
      & 'WSCP=   ',f10.6,' (SC-p)'/
      & 'WELEC=  ',f10.6,' (p-p electr)'/
      & 'WVDWPP= ',f10.6,' (p-p VDW)'/
index ee55c93..00862f2 100644 (file)
@@ -5910,6 +5910,7 @@ C
       logical lprn /.false./, lprn1 /.false./
       etheta=0.0D0
       do i=ithet_start,ithet_end
+        if (i.le.2) cycle
 c        print *,i,itype(i-1),itype(i),itype(i-2)
         if ((itype(i-1).eq.ntyp1).or.itype(i-2).eq.ntyp1
      &  .or.itype(i).eq.ntyp1) cycle
@@ -5926,6 +5927,15 @@ C        print *,i,theta(i)
           sinkt(k)=dsin(k*theti2)
         enddo
 C        print *,ethetai
+        if (i.eq.3) then
+          phii=0.0d0
+          ityp1=nthetyp+1
+          do k=1,nsingle
+            cosph1(k)=0.0d0
+            sinph1(k)=0.0d0
+          enddo
+        else
+
         if (i.gt.3 .and. itype(i-3).ne.ntyp1) then
 #ifdef OSF
           phii=phi(i)
@@ -5947,6 +5957,7 @@ C propagation of chirality for glycine type
             sinph1(k)=0.0d0
           enddo 
         endif
+        endif
         if (i.lt.nres .and. itype(i+1).ne.ntyp1) then
 #ifdef OSF
           phii1=phi(i+1)
index a5d25b3..9347cfa 100644 (file)
@@ -743,11 +743,36 @@ c------------------------------------------------------------------------------
       endif
       call int_from_cart1(.false.)
       do j=nnt+1,nct
+        if (wliptran.gt.0d0) then
+        if (itype(j-1).ne.ntyp1 .and. itype(j).ne.ntyp1 .and.
+     &    (vbld(j).lt.2.0d0 .or. vbld(j).gt.5.3d0)) then
+          if (iprint.gt.0)
+     &    write (iout,*) "Bad CA-CA bond length",j," ",vbld(j),
+     &      " for conformation",ii,wliptran
+          if (iprint.gt.1) then
+            write (iout,*) "The Cartesian geometry is:"
+            write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
+            write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
+            write (iout,*) "The internal geometry is:"
+            write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
+            write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
+            write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
+            write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
+            write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
+            write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
+          endif
+          if (iprint.gt.0) write (iout,*)
+     &      "This conformation WILL NOT be added to the database."
+          conf_check=.false.
+          return
+        endif
+
+        else
         if (itype(j-1).ne.ntyp1 .and. itype(j).ne.ntyp1 .and. 
      &    (vbld(j).lt.2.0d0 .or. vbld(j).gt.5.0d0)) then
           if (iprint.gt.0) 
      &    write (iout,*) "Bad CA-CA bond length",j," ",vbld(j),
-     &      " for conformation",ii
+     &      " for conformation",ii,wliptran
           if (iprint.gt.1) then
             write (iout,*) "The Cartesian geometry is:"
             write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
@@ -765,6 +790,7 @@ c------------------------------------------------------------------------------
           conf_check=.false.
           return
         endif
+        endif
       enddo
       do j=nnt,nct
         itj=itype(j)
index be6ed2a..544a6c7 100644 (file)
@@ -2016,7 +2016,7 @@ C      write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e
         gcorr_loc(i)=0.0d0
       enddo
       do i=iatel_s,iatel_e
-          write (iout,*) i,"i2",itype(i)
+C          write (iout,*) i,"i2",itype(i)
           if (i.eq.1) cycle 
 C           if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1
 C     &  .or. itype(i+2).eq.ntyp1) cycle
@@ -3857,9 +3857,9 @@ C
       etheta=0.0D0
 c      write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1)
       do i=ithet_start,ithet_end
-C         if (i.eq.2) cycle
+         if (i.eq.2) cycle
 C        if (itype(i-1).eq.ntyp1) cycle
-        if (i.le.2) cycle
+C        if (i.le.2) cycle
         if ((itype(i-1).eq.ntyp1).or.itype(i-2).eq.ntyp1
      &  .or.itype(i).eq.ntyp1) cycle
         if (iabs(itype(i+1)).eq.20) iblock=2