update
[unres.git] / source / wham / src-M / readrtns.F
index 36c13b1..1bf78e8 100644 (file)
@@ -17,6 +17,9 @@
       include "COMMON.FREE"
       include "COMMON.CONTROL"
       include "COMMON.ENERGIES"
+      include "COMMON.SPLITELE"
+      include "COMMON.SBRIDGE"
+      include "COMMON.SHIELD"
       character*800 controlcard
       integer i,j,k,ii,n_ene_found
       integer ind,itype1,itype2,itypf,itypsc,itypp
@@ -25,7 +28,7 @@
       character*16 ucase
       character*16 key
       external ucase
-
+      double precision pi
       call card_concat(controlcard,.true.)
       call readi(controlcard,"N_ENE",n_ene,max_ene)
       if (n_ene.gt.max_ene) then
       hamil_rep=index(controlcard,"HAMIL_REP").gt.0
       write (iout,*) "Number of energy parameter sets",nparmset
       call multreadi(controlcard,"ISAMPL",isampl,nparmset,1)
+      do i=1,nparmset
+        if (isampl(i).eq.0) then
+          write (iout,*) "ERROR: isampl is 0 for parmset",i
+          call flush(iout)
+          stop
+        endif
+      enddo
       write (iout,*) "MaxSlice",MaxSlice
       call readi(controlcard,"NSLICE",nslice,1)
       call flush(iout)
@@ -66,6 +76,8 @@
         return1
       endif
       indpdb=0
+      energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
+      rmsrgymap = (index(controlcard,'RMSRGYMAP').gt.0)
       if (index(controlcard,"CLASSIFY").gt.0) indpdb=1
       call reada(controlcard,"DELTA",delta,1.0d-2)
       call readi(controlcard,"EINICHECK",einicheck,2)
       call reada(controlcard,"DELTRGY",deltrgy,5.0d-2)
       call readi(controlcard,"RESCALE",rescale_mode,1)
       check_conf=index(controlcard,"NO_CHECK_CONF").eq.0
+      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,'DISTCHAINMAX',distchainmax,50.0d0)
+       call reada(controlcard,'BOXX',boxxsize,100.0d0)
+       call reada(controlcard,'BOXY',boxysize,100.0d0)
+       call reada(controlcard,'BOXZ',boxzsize,100.0d0)
+c Cutoff range for interactions
+      call reada(controlcard,"R_CUT",r_cut,15.0d0)
+      call reada(controlcard,"LAMBDA",rlamb,0.3d0)
+      call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
+      call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
+      if (lipthick.gt.0.0d0) then
+       bordliptop=(boxzsize+lipthick)/2.0
+       bordlipbot=bordliptop-lipthick
+C      endif
+      if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
+     & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
+      buflipbot=bordlipbot+lipbufthick
+      bufliptop=bordliptop-lipbufthick
+      if ((lipbufthick*2.0d0).gt.lipthick)
+     &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
+      endif
+      write(iout,*) "bordliptop=",bordliptop
+      write(iout,*) "bordlipbot=",bordlipbot
+      write(iout,*) "bufliptop=",bufliptop
+      write(iout,*) "buflipbot=",buflipbot
       call readi(controlcard,'SYM',symetr,1)
       write (iout,*) "DISTCHAINMAX",distchainmax
       write (iout,*) "delta",delta
       zscfile=index(controlcard,"ZSCFILE").gt.0
       with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
       write (iout,*) "with_dihed_constr ",with_dihed_constr
+      with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
+      write (iout,*) "with_theta_constr ",with_theta_constr
+      call readi(controlcard,'SHIELD',shield_mode,0)
+        write(iout,*) "shield_mode",shield_mode
+C      endif
+      call readi(controlcard,'TORMODE',tor_mode,0)
+        write(iout,*) "torsional and valence angle mode",tor_mode
+      if (shield_mode.gt.0) then
+      pi=3.141592d0
+C VSolvSphere the volume of solving sphere
+C      print *,pi,"pi"
+C rpp(1,1) is the energy r0 for peptide group contact and will be used for it 
+C there will be no distinction between proline peptide group and normal peptide
+C group in case of shielding parameters
+      VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
+      VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
+      write (iout,*) VSolvSphere,VSolvSphere_div
+C long axis of side chain 
+C      do i=1,ntyp
+C      long_r_sidechain(i)=vbldsc0(1,i)
+C      short_r_sidechain(i)=sigma0(i)
+C      enddo
+      buff_shield=1.0d0
+      endif
+
       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
+      dyn_ss=index(controlcard,"DYN_SS").gt.0
+      adaptive = index(controlcard,"ADAPTIVE").gt.0
       return
       end
 c------------------------------------------------------------------------------
@@ -400,7 +465,7 @@ c-------------------------------------------------------------------------------
       external ilen,iroof
       double precision rmsdev,energia(0:max_ene),efree,eini,temp
       double precision prop(maxQ)
-      integer ntot_all(maxslice,0:maxprocs-1)
+      integer ntot_all(maxslice,0:maxprocs-1), maxslice_buff
       integer iparm,ib,iib,ir,nprop,nthr,npars
       double precision etot,time
       integer ixdrf,iret 
@@ -531,7 +596,13 @@ c DA scratchfile.
 
 #ifdef MPI
 c Check if everyone has the same number of conformations
-      call MPI_Allgather(stot(1),maxslice,MPI_INTEGER,
+
+c      call MPI_ALLgather(MPI_IN_PLACE,stot(1),MPI_DATATYPE_NULL,
+c     &  ntot_all(1,0),maxslice,MPI_INTEGER,MPI_Comm_World,IERROR)
+
+      maxslice_buff=maxslice
+
+      call MPI_Allgather(stot(1),maxslice_buff,MPI_INTEGER,
      &  ntot_all(1,0),maxslice,MPI_INTEGER,MPI_Comm_World,IERROR)
       lerr=.false.
       do i=0,nprocs-1
@@ -772,3 +843,156 @@ c-------------------------------------------------------------------------------
       iroof = ii
       return
       end
+c-------------------------------------------------------------------------------
+      subroutine read_dist_constr
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.CONTROL'
+      include 'COMMON.CHAIN'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.SBRIDGE'
+      integer ifrag_(2,100),ipair_(2,100)
+      double precision wfrag_(100),wpair_(100)
+      character*500 controlcard
+      logical normalize
+      print *, "WCHODZE"
+      write (iout,*) "Calling read_dist_constr"
+c      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
+c      call flush(iout)
+      call card_concat(controlcard,.true.)
+      call readi(controlcard,"NFRAG",nfrag_,0)
+      call readi(controlcard,"NPAIR",npair_,0)
+      call readi(controlcard,"NDIST",ndist_,0)
+      call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
+      call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
+      call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
+      call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
+      call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
+      normalize = index(controlcard,"NORMALIZE").gt.0
+c      write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
+c      write (iout,*) "IFRAG"
+c      do i=1,nfrag_
+c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
+c      enddo
+c      write (iout,*) "IPAIR"
+c      do i=1,npair_
+c        write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
+c      enddo
+      call flush(iout)
+      do i=1,nfrag_
+        if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
+        if (ifrag_(2,i).gt.nstart_sup+nsup-1)
+     &    ifrag_(2,i)=nstart_sup+nsup-1
+c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
+        call flush(iout)
+        if (wfrag_(i).gt.0.0d0) then
+        do j=ifrag_(1,i),ifrag_(2,i)-1
+          do k=j+1,ifrag_(2,i)
+c            write (iout,*) "j",j," k",k
+            ddjk=dist(j,k)
+            if (constr_dist.eq.1) then
+            nhpb=nhpb+1
+            ihpb(nhpb)=j
+            jhpb(nhpb)=k
+              dhpb(nhpb)=ddjk
+            forcon(nhpb)=wfrag_(i) 
+            else if (constr_dist.eq.2) then
+              if (ddjk.le.dist_cut) then
+                nhpb=nhpb+1
+                ihpb(nhpb)=j
+                jhpb(nhpb)=k
+                dhpb(nhpb)=ddjk
+                forcon(nhpb)=wfrag_(i) 
+              endif
+            else
+              nhpb=nhpb+1
+              ihpb(nhpb)=j
+              jhpb(nhpb)=k
+              dhpb(nhpb)=ddjk
+              forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
+            endif
+            write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
+     &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
+          enddo
+        enddo
+        endif
+      enddo
+      do i=1,npair_
+        if (wpair_(i).gt.0.0d0) then
+        ii = ipair_(1,i)
+        jj = ipair_(2,i)
+        if (ii.gt.jj) then
+          itemp=ii
+          ii=jj
+          jj=itemp
+        endif
+        do j=ifrag_(1,ii),ifrag_(2,ii)
+          do k=ifrag_(1,jj),ifrag_(2,jj)
+            nhpb=nhpb+1
+            ihpb(nhpb)=j
+            jhpb(nhpb)=k
+            forcon(nhpb)=wpair_(i)
+            dhpb(nhpb)=dist(j,k)
+            write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
+     &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
+          enddo
+        enddo
+        endif
+      enddo 
+      write (iout,*) "Distance restraints as read from input"
+      do i=1,ndist_
+        if (constr_dist.eq.11) then
+          read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
+     &     ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
+c        fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
+          if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle
+          nhpb=nhpb+1
+          write (iout,'(a,4i5,2f8.2,2f10.5)') "+dist.restr ",
+     &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
+     &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb)
+          if (ibecarb(i).gt.0) then
+            ihpb(i)=ihpb(i)+nres
+            jhpb(i)=jhpb(i)+nres
+          endif
+        else
+C        print *,"in else"
+          read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
+     &     ibecarb(i),forcon(nhpb+1)
+          if (forcon(nhpb+1).gt.0.0d0) then
+          nhpb=nhpb+1
+          if (ibecarb(i).gt.0) then
+            ihpb(i)=ihpb(i)+nres
+            jhpb(i)=jhpb(i)+nres
+          endif
+          if (dhpb(nhpb).eq.0.0d0)
+     &       dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
+        endif
+          write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
+     &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
+        endif
+C        read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
+C        if (forcon(nhpb+1).gt.0.0d0) then
+C          nhpb=nhpb+1
+C          dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
+      enddo
+      if (constr_dist.eq.11 .and. normalize) then
+        fordepthmax=fordepth(1)
+        do i=2,nhpb
+          if (fordepth(i).gt.fordepthmax) fordepthmax=fordepth(i)
+        enddo
+        do i=1,nhpb
+          fordepth(i)=fordepth(i)/fordepthmax
+        enddo
+        write (iout,'(/a/4a5,2a8,2a10)')
+     &   "Normalized Lorenzian-like distance restraints",
+     &   "   Nr"," res1"," res2"," beta","   d1","   d2","    k","    V"
+        do i=1,nhpb
+          write (iout,'(4i5,2f8.2,2f10.5)')i,ihpb(i),jhpb(i),ibecarb(i),
+     &      dhpb(i),dhpb1(i),forcon(i),fordepth(i)
+        enddo
+      endif
+      write (iout,*) 
+      call hpb_partition
+      call flush(iout)
+      return
+      end