changes in wham and cluser + cutoff corr
[unres4.git] / source / unres / control.F90
index 86f23f9..a1e49d3 100644 (file)
       itordp_nucl= 130
 !      ielep_nucl= 131
       isidep_nucl=132
-
-
+      iscpp_nucl=133
+      isidep_scbase=141
+      isidep_pepbase=142
+      isidep_scpho=143
+      isidep_peppho=144
 
       iliptranpar=60
       itube=61
+!     IONS
+      iion=401
 #if defined(WHAM_RUN) || defined(CLUSTER)
 !
 ! setting the mpi variables for WHAM
 
 !el local variables
       integer :: ind_scint=0,ind_scint_old,ii,jj,i,j,iint,itmp
-
+      integer :: ind_scint_nucl=0
 #ifdef MPI
       integer :: my_sc_int(0:nfgtasks-1),my_ele_int(0:nfgtasks-1)
       integer :: my_sc_intt(0:nfgtasks),my_ele_intt(0:nfgtasks)
       integer :: n_sc_int_tot,my_sc_inde,my_sc_inds,ind_sctint,npept
+      integer :: n_sc_int_tot_nucl,my_sc_inde_nucl,my_sc_inds_nucl, &
+         ind_sctint_nucl,npept_nucl
+
       integer :: nele_int_tot,my_ele_inds,my_ele_inde,ind_eleint_old,&
             ind_eleint,ijunk,nele_int_tot_vdw,my_ele_inds_vdw,&
             my_ele_inde_vdw,ind_eleint_vdw,ind_eleint_vdw_old,&
             ind_scpint_old,nsumgrad,nlen,ngrad_start,ngrad_end,&
             ierror,k,ierr,iaux,ncheck_to,ncheck_from,ind_typ,&
             ichunk,int_index_old
-      integer,dimension(5) :: nct_molec,nnt_molec
+      integer :: nele_int_tot_nucl,my_ele_inds_nucl,my_ele_inde_nucl,&
+            ind_eleint_old_nucl,ind_eleint_nucl,nele_int_tot_vdw_nucl,&
+            my_ele_inds_vdw_nucl,my_ele_inde_vdw_nucl,ind_eleint_vdw_nucl,&
+            ind_eleint_vdw_old_nucl,nscp_int_tot_nucl,my_scp_inds_nucl,&
+            my_scp_inde_nucl,ind_scpint_nucl,ind_scpint_old_nucl,impishi
+!      integer,dimension(5) :: nct_molec,nnt_molec
 !el      allocate(itask_cont_from(0:nfgtasks-1)) !(0:max_fg_procs-1)
 !el      allocate(itask_cont_to(0:nfgtasks-1)) !(0:max_fg_procs-1)
 
       iatsc_e=nct-1
 #endif
       if (iatsc_s.eq.0) iatsc_s=1
+!----------------- scaling for nucleic acid GB
+      n_sc_int_tot_nucl=(nct_molec(2)-nnt_molec(2)+1)*(nct_molec(2)-nnt_molec(2))/2
+      call int_bounds(n_sc_int_tot_nucl,my_sc_inds_nucl,my_sc_inde_nucl)
+!write (iout,*) 'INIT_INT_TABLE nres=',nres,' nnt=',nnt,' nct=',nct
+      if (lprint) &
+        write (iout,*) 'Processor',fg_rank,' CG group',kolor,&
+        ' absolute rank',MyRank,&
+        ' n_sc_int_tot',n_sc_int_tot_nucl,' my_sc_inds=',my_sc_inds_nucl,&
+        ' my_sc_inde',my_sc_inde_nucl
+      ind_sctint_nucl=0
+      iatsc_s_nucl=0
+      iatsc_e_nucl=0
+      do i=1,nres !el   !maxres
+        nint_gr_nucl(i)=0
+        nscp_gr_nucl(i)=0
+        ielstart_nucl(i)=0
+        ielend_nucl(i)=0
+        do j=1,maxint_gr
+          istart_nucl(i,j)=0
+          iend_nucl(i,j)=0
+          iscpstart_nucl(i,j)=0
+          iscpend_nucl(i,j)=0
+        enddo
+      enddo
+      do i=nnt_molec(2),nct_molec(2)-1
+        print*, "inloop2",i
+      call int_partition(ind_scint_nucl,my_sc_inds_nucl,my_sc_inde_nucl,i,&
+           iatsc_s_nucl,iatsc_e_nucl,i+1,nct_molec(2),nint_gr_nucl(i), &
+           istart_nucl(i,1),iend_nucl(i,1),*112)
+        print *,istart_nucl(i,1)
+      enddo
+  112  continue
+       if (iatsc_s_nucl.eq.0) iatsc_s_nucl=1
+       print *,"tu mam",iatsc_s_nucl,iatsc_e_nucl
+
 #ifdef MPI
       if (lprint) write (*,*) 'Processor',fg_rank,' CG Group',kolor,&
          ' absolute rank',myrank,' iatsc_s=',iatsc_s,' iatsc_e=',iatsc_e
         write (iout,'(i3,2(2x,2i3))') &
        i,(istart(i,iint),iend(i,iint),iint=1,nint_gr(i))
       enddo
-      endif
+!      endif
 !      lprint=.false.
+      write (iout,'(a)') 'Interaction array2:' 
+      do i=iatsc_s_nucl,iatsc_e_nucl
+        write (iout,'(i3,2(2x,2i4))') &
+       i,(istart_nucl(i,iint),iend_nucl(i,iint),iint=1,nint_gr_nucl(i))
+      enddo
+      endif
       ispp=4 !?? wham ispp=2
 #ifdef MPI
 ! Now partition the electrostatic-interaction array
-      if (itype(nres_molec(1),1).eq.ntyp1_molec(1)) then
+      if (nres_molec(1).eq.0) then  
+       npept=0
+      elseif (itype(nres_molec(1),1).eq.ntyp1_molec(1)) then
       npept=nres_molec(1)-nnt-1
       else
       npept=nres_molec(1)-nnt
       enddo ! i 
    13 continue
       if (iatel_s.eq.0) iatel_s=1
+!----------now nucleic acid
+!     if (itype(nres_molec(2),2).eq.ntyp1_molec(2)) then
+      npept_nucl=nct_molec(2)-nnt_molec(2)
+!     else
+!     npept_nucl=nct_molec(2)-nnt_molec(2)
+!     endif
+      nele_int_tot_nucl=(npept_nucl-ispp)*(npept_nucl-ispp+1)/2
+      call int_bounds(nele_int_tot_nucl,my_ele_inds_nucl,my_ele_inde_nucl)
+      if (lprint) &
+       write (*,*) 'Processor',fg_rank,' CG group',kolor,&
+        ' absolute rank',MyRank,&
+        ' nele_int_tot',nele_int_tot,' my_ele_inds=',my_ele_inds,&
+                    ' my_ele_inde',my_ele_inde
+      iatel_s_nucl=0
+      iatel_e_nucl=0
+      ind_eleint_nucl=0
+      ind_eleint_old_nucl=0
+!      if (itype(nres_molec(1),1).eq.ntyp1_molec(1)) then
+!      nct_molec(1)=nres_molec(1)-1
+!      else
+!      nct_molec(1)=nres_molec(1)
+!      endif
+!       print *,"nct",nct,nct_molec(1),itype(nres_molec(1),1),ntyp_molec(1)
+      do i=nnt_molec(2),nct_molec(2)-3
+        ijunk=0
+        call int_partition(ind_eleint_nucl,my_ele_inds_nucl,my_ele_inde_nucl,i,&
+          iatel_s_nucl,iatel_e_nucl,i+ispp,nct_molec(2)-1,&
+          ijunk,ielstart_nucl(i),ielend_nucl(i),*113)
+      enddo ! i 
+  113 continue
+      if (iatel_s_nucl.eq.0) iatel_s_nucl=1
+
       nele_int_tot_vdw=(npept-2)*(npept-2+1)/2
 !      write (iout,*) "nele_int_tot_vdw",nele_int_tot_vdw
       call int_bounds(nele_int_tot_vdw,my_ele_inds_vdw,my_ele_inde_vdw)
       enddo ! i 
       if (iatel_s_vdw.eq.0) iatel_s_vdw=1
    15 continue
+      if (iatel_s.eq.0) iatel_s=1
+      if (iatel_s_vdw.eq.0) iatel_s_vdw=1
+      nele_int_tot_vdw_nucl=(npept_nucl-2)*(npept_nucl-2+1)/2
+!      write (iout,*) "nele_int_tot_vdw",nele_int_tot_vdw
+      call int_bounds(nele_int_tot_vdw_nucl,my_ele_inds_vdw_nucl,&
+        my_ele_inde_vdw_nucl)
+!      write (iout,*) "my_ele_inds_vdw",my_ele_inds_vdw,
+!     & " my_ele_inde_vdw",my_ele_inde_vdw
+      ind_eleint_vdw_nucl=0
+      ind_eleint_vdw_old_nucl=0
+      iatel_s_vdw_nucl=0
+      iatel_e_vdw_nucl=0
+      do i=nnt_molec(2),nct_molec(2)-3
+        ijunk=0
+        call int_partition(ind_eleint_vdw_nucl,my_ele_inds_vdw_nucl,&
+          my_ele_inde_vdw_nucl,i,&
+          iatel_s_vdw_nucl,iatel_e_vdw_nucl,i+2,nct_molec(2)-1,&
+          ijunk,ielstart_vdw_nucl(i),&
+          ielend_vdw(i),*115)
+!        write (iout,*) i," ielstart_vdw",ielstart_vdw(i),
+!     &   " ielend_vdw",ielend_vdw(i)
+      enddo ! i 
+      if (iatel_s_vdw.eq.0) iatel_s_vdw_nucl=1
+  115 continue
+
 #else
       iatel_s=nnt
       iatel_e=nct_molec(1)-5 ! ?? wham iatel_e=nct-3
       endif ! lprint
 !     iscp=3
       iscp=2
+      iscp_nucl=2
 ! Partition the SC-p interaction array
 #ifdef MPI
       nscp_int_tot=(npept-iscp+1)*(npept-iscp+1)
         endif
       enddo ! i
    14 continue
+      print *,"before inloop3",iatscp_s,iatscp_e,iscp_nucl
+      nscp_int_tot_nucl=(npept_nucl-iscp_nucl+1)*(npept_nucl-iscp_nucl+1)
+      call int_bounds(nscp_int_tot_nucl,my_scp_inds_nucl,my_scp_inde_nucl)
+      if (lprint) write (iout,*) 'Processor',fg_rank,' CG group',kolor,&
+        ' absolute rank',myrank,&
+        ' nscp_int_tot',nscp_int_tot_nucl,' my_scp_inds=',my_scp_inds_nucl,&
+                    ' my_scp_inde',my_scp_inde_nucl
+      print *,"nscp_int_tot_nucl",nscp_int_tot_nucl,my_scp_inds_nucl,my_scp_inde_nucl
+      iatscp_s_nucl=0
+      iatscp_e_nucl=0
+      ind_scpint_nucl=0
+      ind_scpint_old_nucl=0
+      do i=nnt_molec(2),nct_molec(2)-1
+        print *,"inloop3",i,nnt_molec(2)+iscp,nct_molec(2)-iscp
+        if (i.lt.nnt_molec(2)+iscp) then
+!d        write (iout,*) 'i.le.nnt+iscp'
+          call int_partition(ind_scpint_nucl,my_scp_inds_nucl,&
+            my_scp_inde_nucl,i,iatscp_s_nucl,iatscp_e_nucl,i+iscp,&
+            nct_molec(2),nscp_gr_nucl(i),iscpstart_nucl(i,1),&
+            iscpend_nucl(i,1),*114)
+        else if (i.gt.nct_molec(2)-iscp) then
+!d        write (iout,*) 'i.gt.nct-iscp'
+          call int_partition(ind_scpint_nucl,my_scp_inds_nucl,&
+            my_scp_inde_nucl,i,&
+            iatscp_s_nucl,iatscp_e_nucl,nnt_molec(2),i-iscp,nscp_gr_nucl(i),&
+            iscpstart_nucl(i,1),&
+            iscpend_nucl(i,1),*114)
+        else
+          call int_partition(ind_scpint_nucl,my_scp_inds_nucl,&
+            my_scp_inde_nucl,i,iatscp_s_nucl,iatscp_e_nucl,nnt_molec(2),&
+            i-iscp,nscp_gr_nucl(i),iscpstart_nucl(i,1),&
+           iscpend_nucl(i,1),*114)
+          ii=nscp_gr_nucl(i)+1
+          call int_partition(ind_scpint_nucl,my_scp_inds_nucl,&
+            my_scp_inde_nucl,i,iatscp_s_nucl,iatscp_e_nucl,i+iscp,&
+            nct_molec(2),nscp_gr_nucl(i),iscpstart_nucl(i,ii),&
+            iscpend_nucl(i,ii),*114)
+        endif
+      enddo ! i
+  114 continue
+      print *, "after inloop3",iatscp_s_nucl,iatscp_e_nucl
+      if (iatscp_s_nucl.eq.0) iatscp_s_nucl=1
 #else
       iatscp_s=nnt
       iatscp_e=nct_molec(1)-1
       call int_bounds(nres_molec(1)-2,loc_start,loc_end)
       loc_start=loc_start+1
       loc_end=loc_end+1
+      call int_bounds(nres_molec(2)-2,loc_start_nucl,loc_end_nucl)
+      loc_start_nucl=loc_start_nucl+1+nres_molec(1)
+      loc_end_nucl=loc_end_nucl+1+nres_molec(1)
       call int_bounds(nres_molec(1)-2,ithet_start,ithet_end)
       ithet_start=ithet_start+2
       ithet_end=ithet_end+2
       ibond_nucl_start=ibond_nucl_start+nnt_molec(2)-1
       ibond_nucl_end=ibond_nucl_end+nnt_molec(2)-1
       print *,"NUCLibond",ibond_nucl_start,ibond_nucl_end
+      if (nres_molec(2).ne.0) then
       print *, "before devision",nnt_molec(2),nct_molec(2)-nnt_molec(2)
       call int_bounds(nct_molec(2)-nnt_molec(2),ibondp_nucl_start,ibondp_nucl_end)
       ibondp_nucl_start=ibondp_nucl_start+nnt_molec(2)
       ibondp_nucl_end=ibondp_nucl_end+nnt_molec(2)
+       else
+       ibondp_nucl_start=1
+       ibondp_nucl_end=0
+       endif
       print *,"NUCLibond2",ibondp_nucl_start,ibondp_nucl_end
 
 
        write (*,*) 'Processor:',fg_rank,myrank,' igrad_start',&
          igrad_start,' igrad_end',igrad_end,' ngrad_start',ngrad_start,&
          ' ngrad_end',ngrad_end
-       do i=igrad_start,igrad_end
-         write(*,*) 'Processor:',fg_rank,myrank,i,&
-          jgrad_start(i),jgrad_end(i)
-       enddo
+!       do i=igrad_start,igrad_end
+!         write(*,*) 'Processor:',fg_rank,myrank,i,&
+!          jgrad_start(i),jgrad_end(i)
+!       enddo
       endif
       if (nfgtasks.gt.1) then
         call MPI_Allgather(ivec_start,1,MPI_INTEGER,ivec_displ(0),1,&
         call MPI_Type_contiguous(18,MPI_DOUBLE_PRECISION,MPI_UYZGRAD,&
           IERROR)
         call MPI_Type_commit(MPI_UYZGRAD,IERROR)
+        call MPI_Type_contiguous(maxcontsshi,MPI_INTEGER,MPI_I50,IERROR)
+        call MPI_Type_commit(MPI_I50,IERROR)
+        call MPI_Type_contiguous(maxcontsshi,MPI_DOUBLE_PRECISION,MPI_D50,IERROR)
+        call MPI_Type_commit(MPI_D50,IERROR)
+
+         impishi=maxcontsshi*3
+!        call MPI_Type_contiguous(impishi,MPI_DOUBLE_PRECISION, &
+!        MPI_SHI,IERROR)
+!        call MPI_Type_commit(MPI_SHI,IERROR)
+!        print *,MPI_SHI,"MPI_SHI",MPI_D50
         call MPI_Type_contiguous(2,MPI_DOUBLE_PRECISION,MPI_MU,IERROR)
         call MPI_Type_commit(MPI_MU,IERROR)
         call MPI_Type_contiguous(4,MPI_DOUBLE_PRECISION,MPI_MAT1,IERROR)
 !-----------------------------------------------------------------------------
       subroutine setup_var
 
-      integer :: i
+      integer :: i,mnum
 !      implicit real*8 (a-h,o-z)
 !      include 'DIMENSIONS'
 !      include 'COMMON.IOUNITS'
       nvar=ntheta+nphi
       nside=0
       do i=2,nres-1
+      mnum=molnum(i)
+      write(iout,*) "i",molnum(i)
 #ifdef WHAM_RUN
         if (itype(i,1).ne.10) then
 #else
-        if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
+        if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum) .and. mnum.ne.5) then
 #endif
          nside=nside+1
           ialph(i,1)=nvar+nside
       return
       end subroutine setup_var
 !-----------------------------------------------------------------------------
-! rescode.f
-!-----------------------------------------------------------------------------
-      integer function rescode(iseq,nam,itype,molecule)
-
-      use io_base, only: ucase
-!      implicit real*8 (a-h,o-z)
-!      include 'DIMENSIONS'
-!      include 'COMMON.NAMES'
-!      include 'COMMON.IOUNITS'
-      character(len=3) :: nam  !,ucase
-      integer :: iseq,itype,i
-      integer :: molecule
-      print *,molecule,nam
-      if (molecule.eq.1) then 
-      if (itype.eq.0) then
-
-      do i=-ntyp1_molec(molecule),ntyp1_molec(molecule)
-        if (ucase(nam).eq.restyp(i,molecule)) then
-          rescode=i
-          return
-        endif
-      enddo
-
-      else
-
-      do i=-ntyp1_molec(molecule),ntyp1_molec(molecule)
-        if (nam(1:1).eq.onelet(i)) then
-          rescode=i
-          return  
-        endif  
-      enddo
-
-      endif
-      else if (molecule.eq.2) then
-      do i=1,ntyp1_molec(molecule)
-         print *,nam(1:1),restyp(i,molecule)(1:1) 
-        if (nam(2:2).eq.restyp(i,molecule)(1:1)) then
-          rescode=i
-          return
-        endif
-      enddo
-      else if (molecule.eq.3) then
-       write(iout,*) "SUGAR not yet implemented"
-       stop
-      else if (molecule.eq.4) then
-       write(iout,*) "Explicit LIPID not yet implemented"
-       stop
-      else if (molecule.eq.5) then
-      do i=1,ntyp1_molec(molecule)
-        print *,i,restyp(i,molecule)(1:2)
-        if (ucase(nam(1:2)).eq.restyp(i,molecule)(1:2)) then
-          rescode=i
-          return
-        endif
-      enddo
-      else   
-       write(iout,*) "molecule not defined"
-      endif
-      write (iout,10) iseq,nam
-      stop
-   10 format ('**** Error - residue',i4,' has an unresolved name ',a3)
-      end function rescode
-      integer function sugarcode(sugar,ires)
-      character sugar
-      integer ires
-      if (sugar.eq.'D') then
-        sugarcode=1
-      else if (sugar.eq.' ') then
-        sugarcode=2
-      else
-        write (iout,*) 'UNKNOWN sugar type for residue',ires,' ',sugar
-        stop
-      endif
-      return
-      end function sugarcode
-
-!-----------------------------------------------------------------------------
 ! timing.F
 !-----------------------------------------------------------------------------
 ! $Date: 1994/10/05 16:41:52 $