correcation in ion param and dyn_ss
[unres4.git] / source / cluster / io_clust.F90
index 79cb8d4..ed4ccc3 100644 (file)
       allocate(nss_all(maxstr_proc)) !(maxstr_proc)
       allocate(ihpb_all(maxss,maxstr_proc),jhpb_all(maxss,maxstr_proc))!(maxss,maxstr_proc)
       allocate(iscore(maxconf)) !(maxconf)
-
+      nss_all=0
+      ihpb_all=0
+      jhpb_all=0
 
 #ifdef CHUJ
       ICON=1
 !            call flush(iout)
             if (iret.eq.0) goto 101
             call xdrfint(ixdrf, nss, iret)
-!            write (iout,*) "iret",iret
-!            write (iout,*) "nss",nss
+            write (iout,*) "iret",iret
+!            write (iout,*) "nss",nss,i,"TUTU"
             call flush(iout)
             if (iret.eq.0) goto 101
             do k=1,nss
               if (iret.eq.0) goto 101
               call xdrfint(ixdrf, jhpb(k), iret)
               if (iret.eq.0) goto 101
+!            write(iout,*), ihpb(k),jhpb(k),"TUTU"
             enddo
             call xdrffloat(ixdrf,reini,iret)
             if (iret.eq.0) goto 101
+!            write(iout,*) "TUTU", reini
             call xdrffloat(ixdrf,refree,iret)
+!            write(iout,*) "TUTU", refree
             if (iret.eq.0) goto 101
             call xdrffloat(ixdrf,rmsdev,iret)
             if (iret.eq.0) goto 101
       call int_from_cart1(.false.)
       do j=nnt+1,nct-1
         mnum=molnum(j)
-        write (iout,*) "Check atom",j
-        if (mnum.eq.5) cycle
+!        write (iout,*) "Check atom",j
+        if (mnum.ne.1) cycle
         if (itype(j,mnum).eq.ntyp1_molec(mnum)) cycle
         if (itype(j+1,molnum(j+1)).eq.ntyp1_molec(molnum(j+1))) cycle
 
       do j=nnt,nct
         mnum=molnum(j)
         itj=itype(j,mnum)
-        if (mnum.eq.5) cycle
+        if (mnum.ne.1) cycle
         if (itype(j,1).ne.10 .and. (itype(j,mnum).ne.ntyp1_molec(mnum)) & 
             .and.  (vbld(nres+j)-dsc(iabs(itj))) &
                                   .gt.2.0d0) then
 !------------------------------------------------------------------------------
       subroutine daread_ccoords(istart_conf,iend_conf)
 
+      use energy_data, only: dyn_ss
+
 !      implicit none
 !      include "DIMENSIONS"
 !      include "sizesclu.dat"
         write (iout,*) "Reading binary file, record",iii," ii",ii
         call flush(iout)
 #endif
+        if (dyn_ss) then
+        read(icbase,rec=iii) ((allcart(j,i,ij),j=1,3),i=1,nres), &
+         ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres),&
+         entfac(ii),rmstb(ii)
+        else
         read(icbase,rec=iii) ((allcart(j,i,ij),j=1,3),i=1,nres),&
-          ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres),&
-          nss_all(ij),(ihpb_all(i,ij),jhpb_all(i,ij),i=1,nss),&
-          entfac(ii),rmstb(ii)
+         ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres),&
+         nss_all(ij),(ihpb_all(i,ij),jhpb_all(i,ij),i=1,nss),&
+         entfac(ii),rmstb(ii)
+        endif
+
 #ifdef DEBUG
         write (iout,*) ii,iii,ij,entfac(ii)
         write (iout,'(8f10.5)') ((allcart(j,i,ij),j=1,3),i=1,nres)
 !      implicit none
 !      include "DIMENSIONS"
 !      include "sizesclu.dat"
+      use energy_data, only: dyn_ss
+
 #ifdef MPI
       use MPI_data
       include "mpif.h"
         write (iout,*) "Writing binary file, record",iii," ii",ii
         call flush(iout)
 #endif
+       if (dyn_ss) then
         write(unit_out,rec=iii) ((allcart(j,i,ij),j=1,3),i=1,nres),&
-          ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres),&
-          nss_all(ij),(ihpb_all(i,ij),jhpb_all(i,ij),i=1,nss_all(ij)),&
-          entfac(ii),rmstb(ii)
+         ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres),&
+         entfac(ii),rmstb(ii)
+        else
+        write(unit_out,rec=iii) ((allcart(j,i,ij),j=1,3),i=1,nres), &
+         ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres),&
+         nss_all(ij),(ihpb_all(i,ij),jhpb_all(i,ij),i=1,nss_all(ij)),&
+         entfac(ii),rmstb(ii)
+       endif
+
 #ifdef DEBUG
         write (iout,'(8f10.5)') ((allcart(j,i,ij),j=1,3),i=1,nres)
         write (iout,'(8f10.4)') ((allcart(j,i,ij),j=1,3),i=nnt+nres,&
       call card_concat(controlcard,.true.)
 
       call readi(controlcard,'NRES',nres_molec(1),0)
-      call readi(controlcard,"NRES_NUCL",nres_molec(2),0)
-      call readi(controlcard,"NRES_CAT",nres_molec(5),0)
+      call readi(controlcard,'NRES_NUCL',nres_molec(2),0)
+      call readi(controlcard,'NRES_CAT',nres_molec(5),0)
       nres=0
       do i=1,5
        nres=nres_molec(i)+nres
       enddo
+      print *,"TU",nres_molec(:)
 
 !      call alloc_clust_arrays
       allocate(rcutoff(max_cut+1)) !(max_cut+1)
       call readi(controlcard,'TORMODE',tor_mode,0)
 !C      if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
         write(iout,*) "torsional and valence angle mode",tor_mode
-      call reada(controlcard,"R_CUT_ELE",r_cut_ele,15.0d0)
+      call reada(controlcard,"R_CUT_ELE",r_cut_ele,25.0d0)
       call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
       write(iout,*) "R_CUT_ELE=",r_cut_ele
       call readi(controlcard,'NEND',nend,0)
       character(len=800) :: weightcard
 !      integer rescode
       real(kind=8) :: x(6*nres) !(maxvar)
+      real(kind=8) :: ssscale
       integer  :: itype_pdb(nres) !(maxres)
 !      logical seq_comp
       integer :: i,j,kkk,mnum
       call reada(weightcard,'WCATPROT',wcatprot,0.0d0)
       call reada(weightcard,'TEMP0',temp0,300.0d0)   !!! el
       if (index(weightcard,'SOFT').gt.0) ipot=6
+      call reada(weightcard,'WVDWPP_NUCL',wvdwpp_nucl,0.0D0)
+      call reada(weightcard,'WELPP',welpp,0.0d0)
+      call reada(weightcard,'WVDWPSB',wvdwpsb,0.0d0)
+      call reada(weightcard,'WELPSB',welpsb,0.0D0)
+      call reada(weightcard,'WVDWSB',wvdwsb,0.0d0)
+      call reada(weightcard,'WELSB',welsb,0.0D0)
+      call reada(weightcard,'WBOND_NUCL',wbond_nucl,0.0D0)
+      call reada(weightcard,'WANG_NUCL',wang_nucl,0.0D0)
+      call reada(weightcard,'WSBLOC',wsbloc,0.0D0)
+      call reada(weightcard,'WTOR_NUCL',wtor_nucl,0.0D0)
+!      print *,"KUR..",wtor_nucl
+      call reada(weightcard,'WTORD_NUCL',wtor_d_nucl,0.0D0)
+      call reada(weightcard,'WCORR_NUCL',wcorr_nucl,0.0D0)
+      call reada(weightcard,'WCORR3_NUC',wcorr3_nucl,0.0D0)
 ! 12/1/95 Added weight for the multi-body term WCORR
       call reada(weightcard,'WCORRH',wcorr,1.0D0)
+      call reada(weightcard,'WSCBASE',wscbase,0.0D0)
+      call reada(weightcard,'WPEPBASE',wpepbase,1.0d0)
+      call reada(weightcard,'WSCPHO',wscpho,0.0d0)
+      call reada(weightcard,'WPEPPHO',wpeppho,0.0d0)
+
+      call reada(weightcard,"D0CM",d0cm,3.78d0)
+      call reada(weightcard,"AKCM",akcm,15.1d0)
+      call reada(weightcard,"AKTH",akth,11.0d0)
+      call reada(weightcard,"AKCT",akct,12.0d0)
+      call reada(weightcard,"V1SS",v1ss,-1.08d0)
+      call reada(weightcard,"V2SS",v2ss,7.61d0)
+      call reada(weightcard,"V3SS",v3ss,13.7d0)
+      call reada(weightcard,"EBR",ebr,-5.50D0)
+      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)
+      call reada(weightcard,"SSSCALE",ssscale,1.0D0)
+      dyn_ss=(index(weightcard,'DYN_SS').gt.0)
+
+      call reada(weightcard,"HT",Ht,0.0D0)
+      if (dyn_ss) then
+       ss_depth=(ebr/wsc-0.25*eps(1,1))*ssscale
+        Ht=(Ht/wsc-0.25*eps(1,1))*ssscale
+        akcm=akcm/wsc*ssscale
+        akth=akth/wsc*ssscale
+        akct=akct/wsc*ssscale
+        v1ss=v1ss/wsc*ssscale
+        v2ss=v2ss/wsc*ssscale
+        v3ss=v3ss/wsc*ssscale
+      else
+        ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
+      endif
+
       if (wcorr4.gt.0.0d0) wcorr=wcorr4
       weights(1)=wsc
       weights(2)=wscp
       weights(18)=scal14
 !el      weights(19)=wsccor !!!!!!!!!!!!!!!!
       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
+          weights(41)=wcatcat
+          weights(42)=wcatprot
+          weights(46)=wscbase
+          weights(47)=wscpho
+          weights(48)=wpeppho
+
+
       write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,&
         wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,&
         wturn4,wturn6,wsccor
       else
         read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
       endif
+      print *,nres_molec(:),nres
 ! Convert sequence to numeric code
       do i=1,nres_molec(1)
         mnum=1
       do i=nres_molec(1)+1,nres_molec(1)+nres_molec(2)
         mnum=2
         molnum(i)=2
+        write (iout,*),i,sequence(i)
         itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
+      
       enddo
       do i=nres_molec(1)+nres_molec(2)+1,nres_molec(1)+nres_molec(2)+nres_molec(5)
         mnum=5
@@ -1628,6 +1724,23 @@ write(iout,*)"po ini_int_tab"
 write(iout,*)"przed setup var"
       call setup_var
 write(iout,*)"po setup var"
+      if (ns.gt.0.and.dyn_ss) then
+          do i=nss+1,nhpb
+            ihpb(i-nss)=ihpb(i)
+            jhpb(i-nss)=jhpb(i)
+            forcon(i-nss)=forcon(i)
+            dhpb(i-nss)=dhpb(i)
+          enddo
+          nhpb=nhpb-nss
+          nss=0
+          link_start=1
+          link_end=nhpb
+!          call hpb_partition
+          do i=1,ns
+            dyn_ss_mask(iss(i))=.true.
+          enddo
+      endif
+
       write (iout,*) "molread: REFSTR",refstr
       if (refstr) then
         if (.not.pdbref) then
@@ -1729,6 +1842,8 @@ write(iout,*)"po setup var"
       open (isidep1,file=sidepname,status="old")
       call getenv('SCCORPAR',sccorname)
       open (isccor,file=sccorname,status="old")
+      call getenv('BONDPAR_NUCL',bondname_nucl)
+      open (ibond_nucl,file=bondname_nucl,status='old')
       call getenv('THETPAR_NUCL',thetname_nucl)
       open (ithep_nucl,file=thetname_nucl,status='old')
       call getenv('ROTPAR_NUCL',rotname_nucl)
@@ -1763,6 +1878,9 @@ write(iout,*)"po setup var"
 !
       call getenv('SCPPAR',scpname)
       open (iscpp,file=scpname,status='old')
+      call getenv('SCPPAR_NUCL',scpname_nucl)
+      open (iscpp_nucl,file=scpname_nucl,status='old')
+
 #endif
       return
       end subroutine openunits
@@ -1772,6 +1890,7 @@ write(iout,*)"po setup var"
       subroutine pdboutC(etot,rmsd,tytul)
 
       use energy_data, only: ihpb,jhpb,itype,molnum
+      use energy, only:boxshift,to_box
 !      implicit real*8 (a-h,o-z)
 !      include 'DIMENSIONS'
 !      include 'COMMON.CONTROL'
@@ -1785,22 +1904,30 @@ write(iout,*)"po setup var"
       character(len=50) :: tytul
       character(len=1) :: chainid(10)=(/'A','B','C','D','E','F',&
                                                'G','H','I','J'/)
-      integer :: ica(nres)
+      integer :: ica(nres),k,iti1
       real(kind=8) :: etot,rmsd
       integer :: iatom,ichain,ires,i,j,iti,mnum,mnum1,boxxxshift(3)
-      real(kind=8) :: boxxxx(3)
+      real(kind=8) :: boxxxx(3),cbeg(3),Rdist(3)
       write (ipdb,'(3a,1pe15.5,a,0pf7.2)') 'REMARK ',tytul(:20),&
         ' ENERGY ',etot,' RMS ',rmsd
       iatom=0
       ichain=1
       ires=0
-      boxxxshift(1)=int(c(1,nnt)/boxxsize)
+!      boxxxshift(1)=int(c(1,nnt)/boxxsize)
 !      if (c(1,nnt).lt.0) boxxxshift(1)=boxxxshift(1)-1
-      boxxxshift(2)=int(c(2,nnt)/boxzsize)
+!      boxxxshift(2)=int(c(2,nnt)/boxzsize)
 !      if (c(2,nnt).lt.0) boxxxshift(2)=boxxxshift(2)-1
-      boxxxshift(3)=int(c(3,nnt)/boxzsize)
+!      boxxxshift(3)=int(c(3,nnt)/boxzsize)
 !      if (c(3,nnt).lt.0) boxxxshift(3)=boxxxshift(3)-1
-
+!      do i=1,3
+!        if (boxxxshift(i).gt.2) boxxxshift(i)=2
+!        if (boxxxshift(i).lt.-2) boxxxshift(i)=-2
+!
+!      enddo
+      do j=1,3
+      cbeg(j)=c(j,nnt)
+      enddo
+      call to_box(cbeg(1),cbeg(2),cbeg(3))
       boxxxx(1)=boxxsize
       boxxxx(2)=boxysize
       boxxxx(3)=boxzsize
@@ -1808,25 +1935,37 @@ write(iout,*)"po setup var"
       do i=nnt,nct
         mnum=molnum(i)
         iti=itype(i,mnum)
+        if (i.ne.nct) then
+        iti1=itype(i+1,mnum)
+        mnum1=molnum(i+1)
+        if ((iti.eq.ntyp1_molec(mnum)).and.(iti1.eq.ntyp1_molec(mnum1))) cycle 
+        endif
         if (iti.eq.ntyp1_molec(mnum)) then
           ichain=ichain+1
           ires=0
           write (ipdb,'(a)') 'TER'
         else
         ires=ires+1
+        if (ires.eq.2) then
+         do j=1,3
+         cbeg(j)=c(j,i-1)
+         enddo
+        endif
         iatom=iatom+1
         ica(i)=iatom
-        if (mnum.eq.5) then
-           do j=1,3
-            if ((c(j,i).gt.0).and.(c(j,nnt).lt.0)) then
-           c(j,i)=dmod(c(j,i),boxxxx(j))+(boxxxshift(j)-1)*boxxxx(j)
-            else if ((c(j,i).lt.0).and.(c(j,nnt).gt.0)) then
-           c(j,i)=dmod(c(j,i),boxxxx(j))+(boxxxshift(j)+1)*boxxxx(j)
-             else
-           c(j,i)=dmod(c(j,i),boxxxx(j))+(boxxxshift(j))*boxxxx(j)
-            endif
-           enddo
-        endif
+        call to_box(c(1,i),c(2,i),c(3,i))
+
+        DO k=1,3
+        Rdist(k) = boxshift(c(k,i) - cbeg(k),boxxxx(k))
+        c(k,i)=cbeg(k)+Rdist(k)
+        ENDDO
+      if ((itype(i,molnum(i)).ne.10).and.(molnum(i).ne.5)) then
+        call to_box(c(1,i+nres),c(2,i+nres),c(3,i+nres))
+        DO k=1,3
+        Rdist(k) = boxshift(c(k,i+nres) - cbeg(k),boxxxx(k))
+        c(k,i+nres)=cbeg(k)+Rdist(k)
+        ENDDO
+       endif
         write (ipdb,10) iatom,restyp(iti,mnum),chainid(ichain),&
            ires,(c(j,i),j=1,3),1.0d0,tempfac(1,i)
         if ((iti.ne.10).and.(mnum.ne.5)) then