critical bug fix for ions langvin and fix pdb output for wham and cluster

diff --git a/source/cluster/io_clust.F90 b/source/cluster/io_clust.F90
index ed4ccc3..3b3a134 100644 (file)
--- a/source/cluster/io_clust.F90
+++ b/source/cluster/io_clust.F90
@@ -1493,6 +1493,7 @@
       call reada(weightcard,'WPEPBASE',wpepbase,1.0d0)
       call reada(weightcard,'WSCPHO',wscpho,0.0d0)
       call reada(weightcard,'WPEPPHO',wpeppho,0.0d0)
+      call reada(weightcard,'WCATNUCL',wcatnucl,0.0d0)
 
       call reada(weightcard,"D0CM",d0cm,3.78d0)
       call reada(weightcard,"AKCM",akcm,15.1d0)
@@ -1562,6 +1563,8 @@
           weights(46)=wscbase
           weights(47)=wscpho
           weights(48)=wpeppho
+          weights(49)=wpeppho
+          weights(50)=wcatnucl
 
 
       write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,&
@@ -1870,6 +1873,8 @@ write(iout,*)"po setup var"
       open (itube,file=tubename,status='old')
       call getenv('IONPAR',ionname)
       open (iion,file=ionname,status='old')
+      call getenv('IONPAR_NUCL',ionnuclname)
+      open (iionnucl,file=ionnuclname,status='old')
 
 #ifndef OLDSCP
 !
@@ -1946,7 +1951,7 @@ write(iout,*)"po setup var"
           write (ipdb,'(a)') 'TER'
         else
         ires=ires+1
-        if (ires.eq.2) then
+        if ((ires.eq.2).and.(mnum.ne.5)) then
          do j=1,3
          cbeg(j)=c(j,i-1)
          enddo

diff --git a/source/unres/MD.F90 b/source/unres/MD.F90
index 2b412fd..3bb57e0 100644 (file)
--- a/source/unres/MD.F90
+++ b/source/unres/MD.F90
@@ -2058,6 +2058,7 @@
           do j=1,3    
             adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
             adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
+!            write(iout,*) "adt",adt,"ads",adt2
             dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
             d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
             d_t(j,inres)=d_t_old(j,inres)+adt
@@ -5033,11 +5034,14 @@
 !el       ginvfric(2*nres,2*nres)      !maxres2=2*maxres
 !el      common /przechowalnia/ ginvfric
       
-      logical :: lprn = .false., checkmode = .false.
+      logical :: lprn, checkmode
       integer :: i,j,ind,k,nres2,nres6,mnum
       nres2=2*nres
       nres6=6*nres
-
+      lprn=.false.
+      checkmode=.false.
+!      if (large) lprn=.true.
+!      if (large) checkmode=.true.
       if(.not.allocated(gamvec)) allocate(gamvec(nres6)) !(MAXRES6)
       if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres
       do i=0,nres2

diff --git a/source/unres/MREMD.F90 b/source/unres/MREMD.F90
index 8c4a4b5..003d67b 100644 (file)
--- a/source/unres/MREMD.F90
+++ b/source/unres/MREMD.F90
@@ -411,11 +411,13 @@
            do i=nnt,nct-1
              mnum=(molnum(i))
              stdforcp(i)=stdfp(mnum)*dsqrt(gamp(mnum))
+!             write(iout,*) "stdforcp=",stdforcp(i),itype(i,mnum),i
            enddo
            do i=nnt,nct
             mnum=molnum(i)
              if (itype(i,mnum).ne.ntyp1) stdforcsc(i)=stdfsc(iabs(itype(i,mnum)),mnum)&
                         *dsqrt(gamsc(iabs(itype(i,mnum)),mnum))
+!              write(iout,*) "stdforcsc=",stdforcsc(i),itype(i,mnum),i
            enddo
          endif
        call rescale_weights(t_bath)

diff --git a/source/unres/REMD.F90 b/source/unres/REMD.F90
index fc0c9f1..699631c 100644 (file)
--- a/source/unres/REMD.F90
+++ b/source/unres/REMD.F90
@@ -627,11 +627,11 @@
         ii = ind+m
         mnum=molnum(i)
         iti=itype(i,mnum)
-        if (mnum.eq.5) then
-        mscab=0.0
-        else
+!        if (mnum.eq.5) then
+!        mscab=0.0
+!        else
         mscab=msc(iabs(iti),mnum)
-        endif
+!        endif
         massvec(ii)=mscab
 
         if (iti.ne.10 .and. iti.ne.ntyp1_molec(mnum) .and. mnum.ne.5) then

diff --git a/source/unres/energy.F90 b/source/unres/energy.F90
index d235e01..ec3ffd4 100644 (file)
--- a/source/unres/energy.F90
+++ b/source/unres/energy.F90
@@ -11663,7 +11663,7 @@
         endif
       endif
       endif
-#define DEBUG
+!#define DEBUG
 #ifdef DEBUG
       write (iout,*) "gradc gradx gloc"
       do i=1,nres
@@ -11671,7 +11671,7 @@
          i,(gradc(j,i,icg),j=1,3),(gradx(j,i,icg),j=1,3),gloc(i,icg)
       enddo 
 #endif
-#undef DEBUG
+!#undef DEBUG
 #ifdef TIMING
       time_sumgradient=time_sumgradient+MPI_Wtime()-time01
 #endif
@@ -16730,13 +16730,13 @@ chip1=chip(itypi)
 #endif
 !#define DEBUG
 !el      write (iout,*) "After sum_gradient"
-!#ifdef DEBUG
+#ifdef DEBUG
       write (iout,*) "After sum_gradient"
       do i=1,nres-1
         write (iout,*) i," gradc  ",(gradc(j,i,icg),j=1,3)
         write (iout,*) i," gradx  ",(gradx(j,i,icg),j=1,3)
       enddo
-!#endif
+#endif
 !#undef DEBUG
 ! If performing constraint dynamics, add the gradients of the constraint energy
       if(usampl.and.totT.gt.eq_time) then

diff --git a/source/unres/geometry.F90 b/source/unres/geometry.F90
index c89442f..fe8e89a 100644 (file)
--- a/source/unres/geometry.F90
+++ b/source/unres/geometry.F90
@@ -3398,7 +3398,7 @@
          enddo
       endif 
 !  Settind dE/ddnres-1       
-#define DEBUG
+!#define DEBUG
 #ifdef DEBUG
           j=1
               write(iout,*)"in int to carta",nres-1,gcart(j,nres-1),gloc(nres-3,icg),dphi(j,3,nres), &
@@ -3445,7 +3445,7 @@
               write(iout,*)"in int to cart",i, gxcart(j,i),gloc(ialph(i,1),icg),dalpha(j,3,i), &
               gloc(ialph(i,1)+nside,icg),domega(j,3,i)
 #endif
-#undef DEBUG
+!#undef DEBUG
             enddo
          endif     
        enddo                                                                                                                                                   

diff --git a/source/unres/io_config.F90 b/source/unres/io_config.F90
index 73c18ab..0d92898 100644 (file)
--- a/source/unres/io_config.F90
+++ b/source/unres/io_config.F90
@@ -3293,6 +3293,8 @@
              print *,msc(i,5),restok(i,5)
             enddo
             ip(5)=0.2
+           ! mp(5)=0.2
+           ! pstok(5)=1.0
 !DIR$ NOUNROLL 
       do j=1,ntyp_molec(5)
        do i=1,ntyp
@@ -3935,8 +3937,8 @@
 !          enddo
 !        else
         do j=1,3
-          dcj=(c(j,nres-2)-c(j,nres-3))/2.0
-          c(j,nres)=c(j,nres-1)+dcj
+          dcj=(c(j,nres-nres_molec(5)-2)-c(j,nres-nres_molec(5)-3))/2.0
+          c(j,nres)=c(j,nres-nres_molec(5)-1)+dcj
           c(j,2*nres)=c(j,nres)
         enddo
 !        endif

diff --git a/source/wham/io_wham.F90 b/source/wham/io_wham.F90
index 68457dd..fee1e46 100644 (file)
--- a/source/wham/io_wham.F90
+++ b/source/wham/io_wham.F90
@@ -3908,8 +3908,9 @@ allocate(ww(max_eneW))
 !--------------------------------------------------------------------------------
       subroutine pdboutW(ii,temp,efree,etot,entropy,rmsdev)
 
-      use geometry_data, only:nres,c
+      use geometry_data, only:nres,c,boxxsize,boxysize,boxzsize
       use energy_data, only:nss,nnt,nct,ihpb,jhpb,itype,molnum
+      use energy, only:boxshift
 !      implicit real*8 (a-h,o-z)
 !      include 'DIMENSIONS'
 !      include 'DIMENSIONS.ZSCOPT'
@@ -3924,7 +3925,7 @@ allocate(ww(max_eneW))
                       'D','E','F','G','H','I','J','K','L','M','N','O',&
                       'P','Q','R','S','V','W','X','Y','Z'/),shape(chainid))
       integer,dimension(nres) :: ica !(maxres)
-      real(kind=8) :: temp,efree,etot,entropy,rmsdev
+      real(kind=8) :: temp,efree,etot,entropy,rmsdev,xj,yj,zj
       integer :: ii,i,j,iti,ires,iatom,ichain,mnum
       write(ipdb,'("REMARK CONF",i8," TEMPERATURE",f7.1," RMS",0pf7.2)')&
         ii,temp,rmsdev
@@ -3948,9 +3949,17 @@ allocate(ww(max_eneW))
         ires=ires+1
         iatom=iatom+1
         ica(i)=iatom
+        if (mnum.ne.5) then
         write (ipdb,10) iatom,restyp(iti,mnum),chainid(ichain),&
            ires,(c(j,i),j=1,3)
-        if (iti.ne.10) then
+        else
+        xj=boxshift(c(1,i)-c(1,2),boxxsize)
+        yj=boxshift(c(2,i)-c(2,2),boxysize)
+        zj=boxshift(c(3,i)-c(3,2),boxzsize)
+        write (ipdb,10) iatom,restyp(iti,mnum),chainid(ichain),&
+           ires,c(1,2)+xj,c(2,2)+yj,c(3,2)+zj
+        endif
+        if ((iti.ne.10).and.(mnum.ne.5)) then
           iatom=iatom+1
           write (ipdb,20) iatom,restyp(iti,mnum),chainid(ichain),&
             ires,(c(j,nres+i),j=1,3)