working martini
[unres4.git] / source / unres / control.F90
index 0444b44..599a108 100644 (file)
 !
 ! The following is just to define auxiliary variables used in angle conversion
 !
+!      ifirstrun=0
       pi=4.0D0*datan(1.0D0)
       dwapi=2.0D0*pi
       dwapi3=dwapi/3.0D0
       icsa_in=40
 !rc for ifc error 118
       icsa_pdb=42
+      irotam_end=43
 #endif
       iscpp=25
       icbase=16
       isidep_peppho=144
 
       iliptranpar=60
+      
       itube=61
+!     LIPID MARTINI
+      ilipbond=301
+      ilipnonbond=302
+      imartprot=303 ! this parameters are between protein and martini form of lipid
 !     IONS
       iion=401
+      iionnucl=402
+      iiontran=403 ! this is parameter file for transition metals
+      iwaterwater=404
+      iwatersc=405
 #if defined(WHAM_RUN) || defined(CLUSTER)
 !
 ! setting the mpi variables for WHAM
 !   
 !c     maxfun=5000
 !c     maxit=2000
-      maxfun=500
-      maxit=200
+      maxfun=1000
+      maxit=1000
       tolf=1.0D-2
       rtolf=5.0D-4
 ! 
             nscp_int_tot,my_scp_inds,my_scp_inde,ind_scpint,&
             ind_scpint_old,nsumgrad,nlen,ngrad_start,ngrad_end,&
             ierror,k,ierr,iaux,ncheck_to,ncheck_from,ind_typ,&
-            ichunk,int_index_old
+            ichunk,int_index_old,ibra
       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
+            my_scp_inde_nucl,ind_scpint_nucl,ind_scpint_old_nucl,impishi
+       integer(kind=1),dimension(:,:),allocatable :: remmat
 !      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)
 !... Determine the numbers of start and end SC-SC interaction
 !... to deal with by current processor.
 !write (iout,*) '******INIT_INT_TABLE nres=',nres,' nnt=',nnt,' nct=',nct
+      print *,"in spliting contacts"
       do i=0,nfgtasks-1
         itask_cont_from(i)=fg_rank
         itask_cont_to(i)=fg_rank
       lprint=energy_dec
       itmp=0
       do i=1,5
+       print *,i,nres_molec(i)
        if (nres_molec(i).eq.0) cycle
       itmp=itmp+nres_molec(i)
       if (itype(itmp,i).eq.ntyp1_molec(i)) then
       nnt_molec(i)=itmp+1
       endif
       enddo
+!      if (.not.allocated(nres_molec)) print *,"WHATS WRONG"
       print *,"nres_molec",nres_molec(:)
       print *,"nnt_molec",nnt_molec(:)
       print *,"nct_molec",nct_molec(:)
       iatsc_s=0
       iatsc_e=0
 #endif
+        if(.not.allocated(ielstart_all)) then
 !el       common /przechowalnia/
       allocate(iturn3_start_all(0:nfgtasks))
       allocate(iturn3_end_all(0:nfgtasks))
       allocate(itask_cont_from_all(0:nfgtasks-1,0:nfgtasks-1))
       allocate(itask_cont_to_all(0:nfgtasks-1,0:nfgtasks-1))
 !el----------
+      endif
 !      lprint=.false.
         print *,"NCT",nct_molec(1),nct
       do i=1,nres !el  !maxres
         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
       iturn3_end=iturn3_end+nnt
       iphi_end=iturn3_end+2
       iturn3_start=iturn3_start-1
+      if (iturn3_start.eq.0) iturn3_start=1
       iturn3_end=iturn3_end-1
       call int_bounds(nct_molec(2)-nnt_molec(2)-2,iphi_nucl_start,iphi_nucl_end)
       iphi_nucl_start=iphi_nucl_start+nnt_molec(2)+2
       iphid_end=iturn4_end+2
       iturn4_start=iturn4_start-1
       iturn4_end=iturn4_end-1
+      if (iturn4_start.eq.0) iturn4_start=1
 !      print *,"TUTUTU",nres_molec(1),nres
       call int_bounds(nres_molec(1)-2,ibond_start,ibond_end) 
       ibond_start=ibond_start+1
       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)
+      call int_bounds(nres_molec(2)-1,ibondp_nucl_start,ibondp_nucl_end)
+      ibondp_nucl_start=ibondp_nucl_start+nnt_molec(2)-1
+      ibondp_nucl_end=ibondp_nucl_end+nnt_molec(2)-1
+       else
+       ibondp_nucl_start=1
+       ibondp_nucl_end=0
+       endif
       print *,"NUCLibond2",ibondp_nucl_start,ibondp_nucl_end
 
 
         call int_bounds &
        (ntheta_constr,ithetaconstr_start,ithetaconstr_end)
       endif
+!     HERE MAKING LISTS FOR MARTINI
+      itmp=0
+      do i=1,3
+       itmp=itmp+nres_molec(i)
+      enddo
+!First bonding
+!       call int_bounds(nres_molec(4)-1,ilipbond_start,ilipbond_end)
+       ilipbond_start=1+itmp
+       ilipbond_end=nres_molec(4)-1+itmp
+!angles
+       call int_bounds(nres_molec(4)-1,ilipbond_start_tub,ilipbond_end_tub)
+       ilipbond_start_tub=1+itmp
+       ilipbond_end_tub=nres_molec(4)-1+itmp
 
+!       call int_bounds(nres_molec(4)-2,ilipang_start,ilipang_end)
+       ilipang_start=2+itmp
+       ilipang_end=itmp+nres_molec(4)-1
+!      create LJ LIST MAXIMUM
+!      Eliminate branching from list
+       if(.not.allocated(remmat))&
+        allocate(remmat(itmp+1:nres_molec(4)+itmp,itmp+1:nres_molec(4)+itmp))
+          remmat=0
+       do i=1+itmp,nres_molec(4)-1+itmp
+        if (itype(i,4).eq.12) ibra=i
+        if (itype(i,4).eq.ntyp1_molec(4)-1) then
+!        remmat(ibra-1,i+1)=1
+        remmat(ibra,i+1)=1
+!        remmat(ibra+1,i+1)=1
+        endif
+       enddo
+       maxljliplist=0
+       if (.not.allocated(mlipljlisti)) then
+       allocate (mlipljlisti(nres_molec(4)*nres_molec(4)/2))
+       allocate (mlipljlistj(nres_molec(4)*nres_molec(4)/2))
+       endif
+       do i=1+itmp,nres_molec(4)-1+itmp
+        do j=i+2,nres_molec(4)+itmp
+        if ((itype(i,4).le.ntyp_molec(4)).and.(itype(j,4).le.ntyp_molec(4))&
+        .and.(remmat(i,j).eq.0)) then
+        maxljliplist=maxljliplist+1
+        mlipljlisti(maxljliplist)=i
+        mlipljlistj(maxljliplist)=j
+        if (energy_dec) print *,i,j,remmat(i,j),"lj lip list"
+        endif
+        enddo
+       enddo
+!      split the bound of the list
+       call int_bounds(maxljliplist,iliplj_start,iliplj_end)
+       iliplj_start=iliplj_start
+       iliplj_end=iliplj_end
+!      now the electrostatic list
+       maxelecliplist=0
+       if (.not.allocated(mlipeleclisti)) then
+       allocate (mlipeleclisti(nres_molec(4)*nres_molec(4)/2))
+       allocate (mlipeleclistj(nres_molec(4)*nres_molec(4)/2))
+       endif
+       do i=1+itmp,nres_molec(4)-1+itmp
+        do j=i+2,nres_molec(4)+itmp
+        if ((itype(i,4).le.4).and.(itype(j,4).le.4)) then
+        maxelecliplist=maxelecliplist+1
+        mlipeleclisti(maxelecliplist)=i
+        mlipeleclistj(maxelecliplist)=j
+        endif
+        enddo
+       enddo
+       call int_bounds(maxelecliplist,ilip_elec_start,ilipelec_end)
+       ilip_elec_start=ilip_elec_start
+       ilipelec_end=ilipelec_end
 !      nsumgrad=(nres-nnt)*(nres-nnt+1)/2
 !      nlen=nres-nnt+1
       nsumgrad=(nres-nnt)*(nres-nnt+1)/2
         jgrad_start(i)=i+1
         jgrad_end(i)=nres
       enddo
+! THIS SHOULD BE FIXED
+      itmp=0
+      do i=1,4
+       itmp=itmp+nres_molec(i)
+      enddo
+      call int_bounds(nres_molec(5),icatb_start,icatb_end)
+      icatb_start=icatb_start+itmp
+      icatb_end=icatb_end+itmp
+
+
+
       if (lprint) then 
         write (*,*) 'Processor:',fg_rank,' CG group',kolor,&
        ' absolute rank',myrank,&
         write (iout,*) "iturn4_end_all",&
           (iturn4_end_all(i),i=0,nfgtasks-1)
         write (iout,*) "The ielstart_all array"
+!        do i=0,nfgtasks-1
+!         if (iturn3_start_all(i).le.0) iturn3_start_all(i)=1
+!         if (iturn4_start_all(i).le.0) iturn4_start_all(i)=1
+!        enddo
         do i=nnt,nct
           write (iout,'(20i4)') i,(ielstart_all(i,j),j=0,nfgtasks-1)
         enddo
 !        write (iout,*) "Loop backward iatel_e-1",iatel_e-1,
 !     &     " iatel_e",iatel_e
 !        call flush(iout)
+#ifndef NEWCORR
         nat_sent=0
         do i=iatel_s,iatel_e
 !          write (iout,*) "i",i," ielstart",ielstart(i),
             iat_sent(nat_sent)=i
           endif
         enddo
+#endif
         if (lprint) then
         write (iout,*)"After longrange ntask_cont_from",ntask_cont_from,&
          " ntask_cont_to",ntask_cont_to
         write (iout,*) "itask_cont_to",&
           (itask_cont_to(i),i=1,ntask_cont_to)
         call flush(iout)
+#ifndef NEWCORR
         write (iout,*) "iint_sent"
         do i=1,nat_sent
           ii=iat_sent(i)
           write (iout,'(20i4)') ii,(j,(iint_sent(k,j,i),k=1,4),&
             j=ielstart(ii),ielend(ii))
         enddo
+#endif
         write (iout,*) "iturn3_sent iturn3_start",iturn3_start,&
           " iturn3_end",iturn3_end
         write (iout,'(20i4)') (i,(iturn3_sent(j,i),j=1,4),&
           itask_cont_from(0),CONT_FROM_GROUP,IERR)
         call MPI_Group_incl(fg_group,ntask_cont_to+1,itask_cont_to(0),&
           CONT_TO_GROUP,IERR)
+#ifndef NEWCORR
         do i=1,nat_sent
           ii=iat_sent(i)
           iaux=4*(ielend(ii)-ielstart(ii)+1)
+          if (iaux.lt.0) iaux=0 
           call MPI_Group_translate_ranks(fg_group,iaux,&
             iint_sent(1,ielstart(ii),i),CONT_TO_GROUP,&
             iint_sent_local(1,ielstart(ii),i),IERR )
 !          write (iout,*) "Ranks translated i=",i
 !          call flush(iout)
         enddo
+#endif
         iaux=4*(iturn3_end-iturn3_start+1)
+          if (iaux.lt.0) iaux=0
         call MPI_Group_translate_ranks(fg_group,iaux,&
            iturn3_sent(1,iturn3_start),CONT_TO_GROUP,&
            iturn3_sent_local(1,iturn3_start),IERR)
         iaux=4*(iturn4_end-iturn4_start+1)
+          if (iaux.lt.0) iaux=0
         call MPI_Group_translate_ranks(fg_group,iaux,&
            iturn4_sent(1,iturn4_start),CONT_TO_GROUP,&
            iturn4_sent_local(1,iturn4_start),IERR)
         if (lprint) then
+#ifndef NEWCORR
+
         write (iout,*) "iint_sent_local"
         do i=1,nat_sent
           ii=iat_sent(i)
             j=ielstart(ii),ielend(ii))
           call flush(iout)
         enddo
+#endif
+        if (iturn3_end.gt.0) then
         write (iout,*) "iturn3_sent_local iturn3_start",iturn3_start,&
           " iturn3_end",iturn3_end
         write (iout,'(20i4)') (i,(iturn3_sent_local(j,i),j=1,4),&
            i=iturn4_start,iturn4_end)
         call flush(iout)
         endif
+        endif
         call MPI_Group_free(fg_group,ierr)
         call MPI_Group_free(cont_from_group,ierr)
         call MPI_Group_free(cont_to_group,ierr)
         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)
       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,mnum).ne.ntyp1_molec(mnum) .and. mnum.ne.5) then
+        if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum) .and. mnum.lt.4) 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 $
       end subroutine print_detailed_timing
 #endif
 !-----------------------------------------------------------------------------
+      subroutine homology_partition
+      implicit none
+!      include 'DIMENSIONS'
+!#ifdef MPI
+!      include 'mpif.h'
+!#endif
+!      include 'COMMON.SBRIDGE'
+!      include 'COMMON.IOUNITS'
+!      include 'COMMON.SETUP'
+!      include 'COMMON.CONTROL'
+!      include 'COMMON.INTERACT'
+!      include 'COMMON.HOMOLOGY'
+!d      write(iout,*)"homology_partition: lim_odl=",lim_odl,
+!d     &   " lim_dih",lim_dih
+#ifdef MPI
+      if (me.eq.king .or. .not. out1file) write (iout,*) "MPI"
+      call int_bounds(lim_odl,link_start_homo,link_end_homo)
+      call int_bounds(lim_dih,idihconstr_start_homo, &
+       idihconstr_end_homo)
+      idihconstr_start_homo=idihconstr_start_homo+nnt-1+3
+      idihconstr_end_homo=idihconstr_end_homo+nnt-1+3
+      if (me.eq.king .or. .not. out1file)&
+       write (iout,*) 'Processor',fg_rank,' CG group',kolor,&
+       ' absolute rank',MyRank,&
+       ' lim_odl',lim_odl,' link_start=',link_start_homo,&
+       ' link_end',link_end_homo,' lim_dih',lim_dih,&
+       ' idihconstr_start_homo',idihconstr_start_homo,&
+       ' idihconstr_end_homo',idihconstr_end_homo
+#else
+      write (iout,*) "Not MPI"
+      link_start_homo=1
+      link_end_homo=lim_odl
+      idihconstr_start_homo=nnt+3
+      idihconstr_end_homo=lim_dih+nnt-1+3
+      write (iout,*) &
+       ' lim_odl',lim_odl,' link_start=',link_start_homo, &
+       ' link_end',link_end_homo,' lim_dih',lim_dih,&
+       ' idihconstr_start_homo',idihconstr_start_homo,&
+       ' idihconstr_end_homo',idihconstr_end_homo
+#endif
+      return
+      end subroutine homology_partition
+
 !-----------------------------------------------------------------------------
       end module control