X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fwham%2Fsrc-NEWSC%2Finitialize_p.F;fp=source%2Fwham%2Fsrc-NEWSC%2Finitialize_p.F;h=7ac810959a6032762ff85f35ec5855e43250484a;hb=d101c97dea752458d76055fdbae49c26fff03c1f;hp=0000000000000000000000000000000000000000;hpb=325eda160c9ad2982501e091ca40606a29043712;p=unres.git

diff --git a/source/wham/src-NEWSC/initialize_p.F b/source/wham/src-NEWSC/initialize_p.F
new file mode 100755
index 0000000..7ac8109
--- /dev/null
+++ b/source/wham/src-NEWSC/initialize_p.F
@@ -0,0 +1,577 @@
+      subroutine initialize
+C 
+C Define constants and zero out tables.
+C
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'DIMENSIONS.ZSCOPT'
+#ifdef MPI
+      include 'mpif.h'
+#endif
+      include 'COMMON.IOUNITS'
+      include 'COMMON.CHAIN'
+      include 'COMMON.INTERACT'
+      include 'COMMON.GEO'
+      include 'COMMON.LOCAL'
+      include 'COMMON.TORSION'
+      include 'COMMON.FFIELD'
+      include 'COMMON.SBRIDGE'
+      include 'COMMON.MINIM' 
+      include 'COMMON.DERIV'
+      include "COMMON.WEIGHTS"
+      include "COMMON.NAMES"
+      include "COMMON.TIME1"
+C
+C The following is just to define auxiliary variables used in angle conversion
+C
+      pi=4.0D0*datan(1.0D0)
+      dwapi=2.0D0*pi
+      dwapi3=dwapi/3.0D0
+      pipol=0.5D0*pi
+      deg2rad=pi/180.0D0
+      rad2deg=1.0D0/deg2rad
+      angmin=10.0D0*deg2rad
+C
+C Define I/O units.
+C
+      inp=    1
+      iout=   2
+      ipdbin= 3
+      ipdb=   7
+      imol2=  4
+      igeom=  8
+      intin=  9
+      ithep= 11
+      irotam=12
+      itorp= 13
+      itordp= 23
+      ielep= 14
+      isidep=15 
+      isidep1=22
+      iscpp=25
+      icbase=16
+      ifourier=20
+      istat= 17
+      ientin=18
+      ientout=19
+      ibond=28
+      isccor=29
+C
+C WHAM files
+C
+      ihist=30
+      iweight=31
+      izsc=32
+C
+C Set default weights of the energy terms.
+C
+      wlong=1.0D0
+      welec=1.0D0
+      wtor =1.0D0
+      wang =1.0D0
+      wscloc=1.0D0
+      wstrain=1.0D0
+C
+C Zero out tables.
+C
+      ndih_constr=0
+      do i=1,maxres2
+	do j=1,3
+	  c(j,i)=0.0D0
+	  dc(j,i)=0.0D0
+        enddo
+      enddo
+      do i=1,maxres
+	do j=1,3
+	  xloc(j,i)=0.0D0
+        enddo
+      enddo
+      do i=1,ntyp
+	do j=1,ntyp
+	  aa(i,j)=0.0D0
+	  bb(i,j)=0.0D0
+	  augm(i,j)=0.0D0
+	  sigma(i,j)=0.0D0
+	  r0(i,j)=0.0D0
+	  chi(i,j)=0.0D0
+        enddo
+	do j=1,2
+	  bad(i,j)=0.0D0
+        enddo
+	chip(i)=0.0D0
+	alp(i)=0.0D0
+	sigma0(i)=0.0D0
+	sigii(i)=0.0D0
+	rr0(i)=0.0D0
+	a0thet(i)=0.0D0
+	do j=1,2
+	  athet(j,i)=0.0D0
+	  bthet(j,i)=0.0D0
+        enddo
+	do j=0,3
+	  polthet(j,i)=0.0D0
+        enddo
+	do j=1,3
+	  gthet(j,i)=0.0D0
+        enddo
+	theta0(i)=0.0D0
+	sig0(i)=0.0D0
+	sigc0(i)=0.0D0
+	do j=1,maxlob
+	  bsc(j,i)=0.0D0
+	  do k=1,3
+	    censc(k,j,i)=0.0D0
+          enddo
+          do k=1,3
+	    do l=1,3
+	      gaussc(l,k,j,i)=0.0D0
+            enddo
+          enddo
+	  nlob(i)=0
+        enddo
+      enddo
+      nlob(ntyp1)=0
+      dsc(ntyp1)=0.0D0
+      do i=1,maxtor
+	itortyp(i)=0
+	do j=1,maxtor
+	  do k=1,maxterm
+	    v1(k,j,i)=0.0D0
+	    v2(k,j,i)=0.0D0
+          enddo
+        enddo
+      enddo
+      do i=1,maxres
+	itype(i)=0
+	itel(i)=0
+      enddo
+C Initialize the bridge arrays
+      ns=0
+      nss=0 
+      nhpb=0
+      do i=1,maxss
+	iss(i)=0
+      enddo
+      do i=1,maxdim
+	dhpb(i)=0.0D0
+      enddo
+      do i=1,maxres
+	ihpb(i)=0
+	jhpb(i)=0
+      enddo
+C
+C Initialize timing.
+C
+      call set_timers
+C
+C Initialize variables used in minimization.
+C   
+c     maxfun=5000
+c     maxit=2000
+      maxfun=500
+      maxit=200
+      tolf=1.0D-2
+      rtolf=5.0D-4
+C 
+C Initialize the variables responsible for the mode of gradient storage.
+C
+      nfl=0
+      icg=1
+      do i=1,14
+        do j=1,14
+          if (print_order(i).eq.j) then
+            iw(print_order(i))=j
+            goto 1121
+          endif
+        enddo
+1121    continue
+      enddo
+      calc_grad=.false.
+C Set timers and counters for the respective routines
+      t_func = 0.0d0
+      t_grad = 0.0d0
+      t_fhel = 0.0d0
+      t_fbet = 0.0d0
+      t_ghel = 0.0d0
+      t_gbet = 0.0d0
+      t_viol = 0.0d0
+      t_gviol = 0.0d0
+      n_func = 0
+      n_grad = 0
+      n_fhel = 0
+      n_fbet = 0
+      n_ghel = 0
+      n_gbet = 0
+      n_viol = 0
+      n_gviol = 0
+      n_map = 0
+#ifndef SPLITELE
+      nprint_ene=nprint_ene-1
+#endif
+      return
+      end
+c-------------------------------------------------------------------------
+      block data nazwy
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'DIMENSIONS.ZSCOPT'
+      include 'COMMON.NAMES'
+      include 'COMMON.WEIGHTS'
+      include 'COMMON.FFIELD'
+      include 'COMMON.INTERACT'
+      data restyp /
+     &'CYS','MET','PHE','ILE','LEU','VAL','TRP','TYR','ALA','GLY','THR',
+     &'SER','GLN','ASN','GLU','ASP','HIS','ARG','LYS','PRO','D'/
+      data onelet /
+     &'C','M','F','I','L','V','W','Y','A','G','T',
+     &'S','Q','N','E','D','H','R','K','P','X'/
+      data potname /'LJ','LJK','BP','GB','GBV','MM'/
+      data ename / 
+     &   "EVDW SC-SC","EVDW2 SC-p","EES p-p","ECORR4 ","ECORR5 ",
+     &   "ECORR6 ","EELLO ","ETURN3 ","ETURN4 ","ETURN6 ",
+     &   "EBE bend","ESC SCloc","ETORS ","ETORSD ","EHPB","EVDWPP",
+     &   "EVDW2_14","ESTR","ESCCOR","EDIHC","EVDW_T"/
+      data wname /
+     &   "WSC","WSCP","WELEC","WCORR","WCORR5","WCORR6","WEL_LOC",
+     &   "WTURN3","WTURN4","WTURN6","WANG","WSCLOC","WTOR","WTORD",
+     &   "WHPB","WVDWPP","WSCP14","WBOND","WSCCOR","WDIHC","WSC"/
+      data ww0 /1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,
+     &    1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,0.4d0,1.0d0,1.0d0,
+     &    0.0d0,0.0/
+      data nprint_ene /21/
+      data print_order /1,2,3,18,11,12,13,14,4,5,6,7,8,9,10,19,
+     &  16,15,17,20,21/
+c Dielectric constant of water
+      data eps_out /80.0d0/
+      end 
+c---------------------------------------------------------------------------
+      subroutine init_int_table
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'DIMENSIONS.ZSCOPT'
+#ifdef MPI
+      include 'mpif.h'
+#endif
+#ifdef MP
+      include 'COMMON.INFO'
+#endif
+      include 'COMMON.CHAIN'
+      include 'COMMON.INTERACT'
+      include 'COMMON.LOCAL'
+      include 'COMMON.SBRIDGE'
+      include 'COMMON.IOUNITS'
+      logical scheck,lprint
+#ifdef MPL
+      integer my_sc_int(0:max_fg_Procs-1),my_sc_intt(0:max_fg_Procs),
+     & my_ele_int(0:max_fg_Procs-1),my_ele_intt(0:max_fg_Procs)
+C... Determine the numbers of start and end SC-SC interaction 
+C... to deal with by current processor.
+      lprint=.false.
+      if (lprint)
+     &write (iout,*) 'INIT_INT_TABLE nres=',nres,' nnt=',nnt,' nct=',nct
+      n_sc_int_tot=(nct-nnt+1)*(nct-nnt)/2-nss
+      MyRank=MyID-(MyGroup-1)*fgProcs
+      call int_bounds(n_sc_int_tot,my_sc_inds,my_sc_inde)
+      if (lprint)
+     &  write (iout,*) 'Processor',MyID,' MyRank',MyRank,
+     &  ' n_sc_int_tot',n_sc_int_tot,' my_sc_inds=',my_sc_inds,
+     &  ' my_sc_inde',my_sc_inde
+      ind_sctint=0
+      iatsc_s=0
+      iatsc_e=0
+#endif
+      lprint=.false.
+      do i=1,maxres
+        nint_gr(i)=0
+        nscp_gr(i)=0
+        do j=1,maxint_gr
+          istart(i,1)=0
+          iend(i,1)=0
+          ielstart(i)=0
+          ielend(i)=0
+          iscpstart(i,1)=0
+          iscpend(i,1)=0    
+        enddo
+      enddo
+      ind_scint=0
+      ind_scint_old=0
+cd    write (iout,*) 'ns=',ns,' nss=',nss,' ihpb,jhpb',
+cd   &   (ihpb(i),jhpb(i),i=1,nss)
+      do i=nnt,nct-1
+        scheck=.false.
+        do ii=1,nss
+          if (ihpb(ii).eq.i+nres) then
+            scheck=.true.
+            jj=jhpb(ii)-nres
+            goto 10
+          endif
+        enddo
+   10   continue
+cd      write (iout,*) 'i=',i,' scheck=',scheck,' jj=',jj
+        if (scheck) then
+          if (jj.eq.i+1) then
+#ifdef MPL
+            write (iout,*) 'jj=i+1'
+            call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
+     & iatsc_s,iatsc_e,i+2,nct,nint_gr(i),istart(i,1),iend(i,1),*12)
+#else
+            nint_gr(i)=1
+            istart(i,1)=i+2
+            iend(i,1)=nct
+#endif
+          else if (jj.eq.nct) then
+#ifdef MPL
+            write (iout,*) 'jj=nct'
+            call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
+     &  iatsc_s,iatsc_e,i+1,nct-1,nint_gr(i),istart(i,1),iend(i,1),*12)
+#else
+            nint_gr(i)=1
+            istart(i,1)=i+1
+            iend(i,1)=nct-1
+#endif
+          else
+#ifdef MPL
+            call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
+     & iatsc_s,iatsc_e,i+1,jj-1,nint_gr(i),istart(i,1),iend(i,1),*12)
+            ii=nint_gr(i)+1
+            call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
+     & iatsc_s,iatsc_e,jj+1,nct,nint_gr(i),istart(i,ii),iend(i,ii),*12)
+#else
+            nint_gr(i)=2
+            istart(i,1)=i+1
+            iend(i,1)=jj-1
+            istart(i,2)=jj+1
+            iend(i,2)=nct
+#endif
+          endif
+        else
+#ifdef MPL
+          call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
+     &    iatsc_s,iatsc_e,i+1,nct,nint_gr(i),istart(i,1),iend(i,1),*12)
+#else
+          nint_gr(i)=1
+          istart(i,1)=i+1
+          iend(i,1)=nct
+          ind_scint=int_scint+nct-i
+#endif
+        endif
+#ifdef MPL
+        ind_scint_old=ind_scint
+#endif
+      enddo
+   12 continue
+#ifndef MPL
+      iatsc_s=nnt
+      iatsc_e=nct-1
+#endif
+#ifdef MPL
+      if (lprint) then
+        write (iout,*) 'Processor',MyID,' Group',MyGroup
+        write (iout,*) 'iatsc_s=',iatsc_s,' iatsc_e=',iatsc_e
+      endif
+#endif
+      if (lprint) then
+      write (iout,'(a)') 'Interaction array:'
+      do i=iatsc_s,iatsc_e
+        write (iout,'(i3,2(2x,2i3))') 
+     & i,(istart(i,iint),iend(i,iint),iint=1,nint_gr(i))
+      enddo
+      endif
+      ispp=2
+#ifdef MPL
+C Now partition the electrostatic-interaction array
+      npept=nct-nnt
+      nele_int_tot=(npept-ispp)*(npept-ispp+1)/2
+      call int_bounds(nele_int_tot,my_ele_inds,my_ele_inde)
+      if (lprint)
+     & write (iout,*) 'Processor',MyID,' MyRank',MyRank,
+     &  ' nele_int_tot',nele_int_tot,' my_ele_inds=',my_ele_inds,
+     &               ' my_ele_inde',my_ele_inde
+      iatel_s=0
+      iatel_e=0
+      ind_eleint=0
+      ind_eleint_old=0
+      do i=nnt,nct-3
+        ijunk=0
+        call int_partition(ind_eleint,my_ele_inds,my_ele_inde,i,
+     &    iatel_s,iatel_e,i+ispp,nct-1,ijunk,ielstart(i),ielend(i),*13)
+      enddo ! i 
+   13 continue
+#else
+      iatel_s=nnt
+      iatel_e=nct-3
+      do i=iatel_s,iatel_e
+        ielstart(i)=i+2
+        ielend(i)=nct-1
+      enddo
+#endif
+      if (lprint) then
+        write (iout,'(a)') 'Electrostatic interaction array:'
+        do i=iatel_s,iatel_e
+          write (iout,'(i3,2(2x,2i3))') i,ielstart(i),ielend(i)
+        enddo
+      endif ! lprint
+c     iscp=3
+      iscp=2
+C Partition the SC-p interaction array
+#ifdef MPL
+      nscp_int_tot=(npept-iscp+1)*(npept-iscp+1)
+      call int_bounds(nscp_int_tot,my_scp_inds,my_scp_inde)
+      if (lprint)
+     & write (iout,*) 'Processor',MyID,' MyRank',MyRank,
+     &  ' nscp_int_tot',nscp_int_tot,' my_scp_inds=',my_scp_inds,
+     &               ' my_scp_inde',my_scp_inde
+      iatscp_s=0
+      iatscp_e=0
+      ind_scpint=0
+      ind_scpint_old=0
+      do i=nnt,nct-1
+        if (i.lt.nnt+iscp) then
+cd        write (iout,*) 'i.le.nnt+iscp'
+          call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
+     &      iatscp_s,iatscp_e,i+iscp,nct,nscp_gr(i),iscpstart(i,1),
+     &      iscpend(i,1),*14)
+        else if (i.gt.nct-iscp) then
+cd        write (iout,*) 'i.gt.nct-iscp'
+          call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
+     &      iatscp_s,iatscp_e,nnt,i-iscp,nscp_gr(i),iscpstart(i,1),
+     &      iscpend(i,1),*14)
+        else
+          call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
+     &      iatscp_s,iatscp_e,nnt,i-iscp,nscp_gr(i),iscpstart(i,1),
+     &      iscpend(i,1),*14)
+          ii=nscp_gr(i)+1
+          call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
+     &      iatscp_s,iatscp_e,i+iscp,nct,nscp_gr(i),iscpstart(i,ii),
+     &      iscpend(i,ii),*14)
+        endif
+      enddo ! i
+   14 continue
+#else
+      iatscp_s=nnt
+      iatscp_e=nct-1
+      do i=nnt,nct-1
+        if (i.lt.nnt+iscp) then
+          nscp_gr(i)=1
+          iscpstart(i,1)=i+iscp
+          iscpend(i,1)=nct
+        elseif (i.gt.nct-iscp) then
+          nscp_gr(i)=1
+          iscpstart(i,1)=nnt
+          iscpend(i,1)=i-iscp
+        else
+          nscp_gr(i)=2
+          iscpstart(i,1)=nnt
+          iscpend(i,1)=i-iscp
+          iscpstart(i,2)=i+iscp
+          iscpend(i,2)=nct
+        endif 
+      enddo ! i
+#endif
+      if (lprint) then
+        write (iout,'(a)') 'SC-p interaction array:'
+        do i=iatscp_s,iatscp_e
+          write (iout,'(i3,2(2x,2i3))') 
+     &         i,(iscpstart(i,j),iscpend(i,j),j=1,nscp_gr(i))
+        enddo
+      endif ! lprint
+C Partition local interactions
+#ifdef MPL
+      call int_bounds(nres-2,loc_start,loc_end)
+      loc_start=loc_start+1
+      loc_end=loc_end+1
+      call int_bounds(nres-2,ithet_start,ithet_end)
+      ithet_start=ithet_start+2
+      ithet_end=ithet_end+2
+      call int_bounds(nct-nnt-2,iphi_start,iphi_end) 
+      iphi_start=iphi_start+nnt+2
+      iphi_end=iphi_end+nnt+2
+      call int_bounds(nres-3,itau_start,itau_end)
+      itau_start=itau_start+3
+      itau_end=itau_end+3
+      if (lprint) then 
+        write (iout,*) 'Processor:',MyID,
+     & ' loc_start',loc_start,' loc_end',loc_end,
+     & ' ithet_start',ithet_start,' ithet_end',ithet_end,
+     & ' iphi_start',iphi_start,' iphi_end',iphi_end
+        write (*,*) 'Processor:',MyID,
+     & ' loc_start',loc_start,' loc_end',loc_end,
+     & ' ithet_start',ithet_start,' ithet_end',ithet_end,
+     & ' iphi_start',iphi_start,' iphi_end',iphi_end
+      endif
+      if (fgprocs.gt.1 .and. MyID.eq.BossID) then
+        write(iout,'(i10,a,i10,a,i10,a/a,i3,a)') n_sc_int_tot,' SC-SC ',
+     & nele_int_tot,' electrostatic and ',nscp_int_tot,
+     & ' SC-p interactions','were distributed among',fgprocs,
+     & ' fine-grain processors.'
+      endif
+#else
+      loc_start=2
+      loc_end=nres-1
+      ithet_start=3 
+      ithet_end=nres
+      iphi_start=nnt+3
+      iphi_end=nct
+      itau_start=4
+      itau_end=nres
+#endif
+      return
+      end 
+c---------------------------------------------------------------------------
+      subroutine int_partition(int_index,lower_index,upper_index,atom,
+     & at_start,at_end,first_atom,last_atom,int_gr,jat_start,jat_end,*)
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.IOUNITS'
+      integer int_index,lower_index,upper_index,atom,at_start,at_end,
+     & first_atom,last_atom,int_gr,jat_start,jat_end
+      logical lprn
+      lprn=.false.
+      if (lprn) write (iout,*) 'int_index=',int_index
+      int_index_old=int_index
+      int_index=int_index+last_atom-first_atom+1
+      if (lprn) 
+     &   write (iout,*) 'int_index=',int_index,
+     &               ' int_index_old',int_index_old,
+     &               ' lower_index=',lower_index,
+     &               ' upper_index=',upper_index,
+     &               ' atom=',atom,' first_atom=',first_atom,
+     &               ' last_atom=',last_atom
+      if (int_index.ge.lower_index) then
+        int_gr=int_gr+1
+        if (at_start.eq.0) then
+          at_start=atom
+          jat_start=first_atom-1+lower_index-int_index_old
+        else
+          jat_start=first_atom
+        endif
+        if (lprn) write (iout,*) 'jat_start',jat_start
+        if (int_index.ge.upper_index) then
+          at_end=atom
+          jat_end=first_atom-1+upper_index-int_index_old
+          return1
+        else
+          jat_end=last_atom
+        endif
+        if (lprn) write (iout,*) 'jat_end',jat_end
+      endif
+      return
+      end
+c------------------------------------------------------------------------------
+      subroutine hpb_partition
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.SBRIDGE'
+      include 'COMMON.IOUNITS'
+#ifdef MPL
+      include 'COMMON.INFO'
+      call int_bounds(nhpb,link_start,link_end)
+#else
+      link_start=1
+      link_end=nhpb
+#endif
+cd    write (iout,*) 'Processor',MyID,' MyRank',MyRank,
+cd   &  ' nhpb',nhpb,' link_start=',link_start,
+cd   &  ' link_end',link_end
+      return
+      end