source/unres/src_MD-M/initialize_p.F

   1       block data
   2       implicit real*8 (a-h,o-z)
   3       include 'DIMENSIONS'
   4       include 'COMMON.MCM'
   5       include 'COMMON.MD'
   6       data MovTypID
   7      &  /'pool','chain regrow','multi-bond','phi','theta','side chain',
   8      &   'total'/
   9 c Conversion from poises to molecular unit and the gas constant
  10       data cPoise /2.9361d0/, Rb /0.001986d0/
  11       end
  12 c--------------------------------------------------------------------------
  13       subroutine initialize
  14 C
  15 C Define constants and zero out tables.
  16 C
  17       implicit real*8 (a-h,o-z)
  18       include 'DIMENSIONS'
  19 #ifdef MPI
  20       include 'mpif.h'
  21 #endif
  22 #ifndef ISNAN
  23       external proc_proc
  24 #ifdef WINPGI
  25 cMS$ATTRIBUTES C ::  proc_proc
  26 #endif
  27 #endif
  28       include 'COMMON.IOUNITS'
  29       include 'COMMON.CHAIN'
  30       include 'COMMON.INTERACT'
  31       include 'COMMON.GEO'
  32       include 'COMMON.LOCAL'
  33       include 'COMMON.TORSION'
  34       include 'COMMON.FFIELD'
  35       include 'COMMON.SBRIDGE'
  36       include 'COMMON.MCM'
  37       include 'COMMON.MINIM'
  38       include 'COMMON.DERIV'
  39       include 'COMMON.SPLITELE'
  40 c Common blocks from the diagonalization routines
  41       COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(400)
  42       COMMON /MACHSW/ KDIAG,ICORFL,IXDR
  43       logical mask_r
  44 c      real*8 text1 /'initial_i'/
  45
  46       mask_r=.false.
  47 #ifndef ISNAN
  48 c NaNQ initialization
  49       i=-1
  50       rr=dacos(100.0d0)
  51 #ifdef WINPGI
  52       idumm=proc_proc(rr,i)
  53 #else
  54       call proc_proc(rr,i)
  55 #endif
  56 #endif
  57
  58       kdiag=0
  59       icorfl=0
  60       iw=2
  61 C
  62 C The following is just to define auxiliary variables used in angle conversion
  63 C
  64       pi=4.0D0*datan(1.0D0)
  65       dwapi=2.0D0*pi
  66       dwapi3=dwapi/3.0D0
  67       pipol=0.5D0*pi
  68       deg2rad=pi/180.0D0
  69       rad2deg=1.0D0/deg2rad
  70       angmin=10.0D0*deg2rad
  71 C
  72 C Define I/O units.
  73 C
  74       inp=    1
  75       iout=   2
  76       ipdbin= 3
  77       ipdb=   7
  78       icart = 30
  79       imol2=  4
  80       igeom=  8
  81       intin=  9
  82       ithep= 11
  83       ithep_pdb=51
  84       irotam=12
  85       irotam_pdb=52
  86       itorp= 13
  87       itordp= 23
  88       ielep= 14
  89       isidep=15
  90       iscpp=25
  91       icbase=16
  92       ifourier=20
  93       istat= 17
  94       irest1=55
  95       irest2=56
  96       iifrag=57
  97       ientin=18
  98       ientout=19
  99       ibond = 28
 100       isccor = 29
 101 crc for write_rmsbank1
 102       izs1=21
 103 cdr  include secondary structure prediction bias
 104       isecpred=27
 105 C
 106 C CSA I/O units (separated from others especially for Jooyoung)
 107 C
 108       icsa_rbank=30
 109       icsa_seed=31
 110       icsa_history=32
 111       icsa_bank=33
 112       icsa_bank1=34
 113       icsa_alpha=35
 114       icsa_alpha1=36
 115       icsa_bankt=37
 116       icsa_int=39
 117       icsa_bank_reminimized=38
 118       icsa_native_int=41
 119       icsa_in=40
 120 crc for ifc error 118
 121       icsa_pdb=42
 122 C
 123 C Set default weights of the energy terms.
 124 C
 125       wlong=1.0D0
 126       welec=1.0D0
 127       wtor =1.0D0
 128       wang =1.0D0
 129       wscloc=1.0D0
 130       wstrain=1.0D0
 131 C
 132 C Zero out tables.
 133 C
 134 c      print '(a,$)','Inside initialize'
 135 c      call memmon_print_usage()
 136       do i=1,maxres2
 137         do j=1,3
 138           c(j,i)=0.0D0
 139           dc(j,i)=0.0D0
 140         enddo
 141       enddo
 142       do i=1,maxres
 143         do j=1,3
 144           xloc(j,i)=0.0D0
 145         enddo
 146       enddo
 147       do i=1,ntyp
 148         do j=1,ntyp
 149           aa(i,j)=0.0D0
 150           bb(i,j)=0.0D0
 151           augm(i,j)=0.0D0
 152           sigma(i,j)=0.0D0
 153           r0(i,j)=0.0D0
 154           chi(i,j)=0.0D0
 155         enddo
 156         do j=1,2
 157           bad(i,j)=0.0D0
 158         enddo
 159         chip(i)=0.0D0
 160         alp(i)=0.0D0
 161         sigma0(i)=0.0D0
 162         sigii(i)=0.0D0
 163         rr0(i)=0.0D0
 164         a0thet(i)=0.0D0
 165         do j=1,2
 166          do ichir1=-1,1
 167           do ichir2=-1,1
 168           athet(j,i,ichir1,ichir2)=0.0D0
 169           bthet(j,i,ichir1,ichir2)=0.0D0
 170           enddo
 171          enddo
 172         enddo
 173         do j=0,3
 174           polthet(j,i)=0.0D0
 175         enddo
 176         do j=1,3
 177           gthet(j,i)=0.0D0
 178         enddo
 179         theta0(i)=0.0D0
 180         sig0(i)=0.0D0
 181         sigc0(i)=0.0D0
 182         do j=1,maxlob
 183           bsc(j,i)=0.0D0
 184           do k=1,3
 185             censc(k,j,i)=0.0D0
 186           enddo
 187           do k=1,3
 188             do l=1,3
 189               gaussc(l,k,j,i)=0.0D0
 190             enddo
 191           enddo
 192           nlob(i)=0
 193         enddo
 194       enddo
 195       nlob(ntyp1)=0
 196       dsc(ntyp1)=0.0D0
 197       do i=-maxtor,maxtor
 198         itortyp(i)=0
 199 cc      write (iout,*) "TU DOCHODZE",i,itortyp(i)
 200        do iblock=1,2
 201         do j=-maxtor,maxtor
 202           do k=1,maxterm
 203             v1(k,j,i,iblock)=0.0D0
 204             v2(k,j,i,iblock)=0.0D0
 205           enddo
 206         enddo
 207         enddo
 208       enddo
 209       do iblock=1,2
 210        do i=-maxtor,maxtor
 211         do j=-maxtor,maxtor
 212          do k=-maxtor,maxtor
 213           do l=1,maxtermd_1
 214             v1c(1,l,i,j,k,iblock)=0.0D0
 215             v1s(1,l,i,j,k,iblock)=0.0D0
 216             v1c(2,l,i,j,k,iblock)=0.0D0
 217             v1s(2,l,i,j,k,iblock)=0.0D0
 218           enddo !l
 219           do l=1,maxtermd_2
 220            do m=1,maxtermd_2
 221             v2c(m,l,i,j,k,iblock)=0.0D0
 222             v2s(m,l,i,j,k,iblock)=0.0D0
 223            enddo !m
 224           enddo !l
 225         enddo !k
 226        enddo !j
 227       enddo !i
 228       enddo !iblock
 229
 230       do i=1,maxres
 231         itype(i)=0
 232         itel(i)=0
 233       enddo
 234 C Initialize the bridge arrays
 235       ns=0
 236       nss=0
 237       nhpb=0
 238       do i=1,maxss
 239         iss(i)=0
 240       enddo
 241       do i=1,maxdim
 242         dhpb(i)=0.0D0
 243       enddo
 244       do i=1,maxres
 245         ihpb(i)=0
 246         jhpb(i)=0
 247       enddo
 248 C
 249 C Initialize timing.
 250 C
 251       call set_timers
 252 C
 253 C Initialize variables used in minimization.
 254 C
 255 c     maxfun=5000
 256 c     maxit=2000
 257       maxfun=500
 258       maxit=200
 259       tolf=1.0D-2
 260       rtolf=5.0D-4
 261 C
 262 C Initialize the variables responsible for the mode of gradient storage.
 263 C
 264       nfl=0
 265       icg=1
 266 C
 267 C Initialize constants used to split the energy into long- and short-range
 268 C components
 269 C
 270       r_cut=2.0d0
 271       rlamb=0.3d0
 272 #ifndef SPLITELE
 273       nprint_ene=nprint_ene-1
 274 #endif
 275       return
 276       end
 277 c-------------------------------------------------------------------------
 278       block data nazwy
 279       implicit real*8 (a-h,o-z)
 280       include 'DIMENSIONS'
 281       include 'COMMON.NAMES'
 282       include 'COMMON.FFIELD'
 283       data restyp /
 284      &'DD','DAU','DAI','DDB','DSM','DPR','DLY','DAR','DHI','DAS','DGL',
 285      & 'DSG','DGN','DSN','DTH',
 286      &'DYY','DAL','DTY','DTR','DVA','DLE','DIL','DPN','MED','DCY','ZER',
 287      &'CYS','MET','PHE','ILE','LEU','VAL','TRP','TYR','ALA','GLY','THR',
 288      &'SER','GLN','ASN','GLU','ASP','HIS','ARG','LYS','PRO','SME','DBZ',
 289      &'AIB','ABU','D'/
 290       data onelet /
 291      &'z','z','z','z','z','p','k','r','h','d','e','n','q','s','t','g',
 292      &'a','y','w','v','l','i','f','m','c','x',
 293      &'C','M','F','I','L','V','W','Y','A','G','T',
 294      &'S','Q','N','E','D','H','R','K','P','z','z','z','z','X'/
 295       data potname /'LJ','LJK','BP','GB','GBV'/
 296       data ename /
 297      &   "EVDW SC-SC","EVDW2 SC-p","EES p-p","ECORR4 ","ECORR5 ",
 298      &   "ECORR6 ","EELLO ","ETURN3 ","ETURN4 ","ETURN6 ",
 299      &   "EBE bend","ESC SCloc","ETORS ","ETORSD ","EHPB ","EVDWPP ",
 300      &   "ESTR ","EVDW2_14 ","UCONST ", "      ","ESCCOR"/
 301       data wname /
 302      &   "WSC","WSCP","WELEC","WCORR","WCORR5","WCORR6","WEL_LOC",
 303      &   "WTURN3","WTURN4","WTURN6","WANG","WSCLOC","WTOR","WTORD",
 304      &   "WSTRAIN","WVDWPP","WBOND","SCAL14","     ","    ","WSCCOR"/
 305       data nprint_ene /20/
 306       data print_order/1,2,3,11,12,13,14,4,5,6,7,8,9,10,19,18,15,17,16,
 307      & 21,0/
 308       end
 309 c---------------------------------------------------------------------------
 310       subroutine init_int_table
 311       implicit real*8 (a-h,o-z)
 312       include 'DIMENSIONS'
 313 #ifdef MPI
 314       include 'mpif.h'
 315       integer blocklengths(15),displs(15)
 316 #endif
 317       include 'COMMON.CONTROL'
 318       include 'COMMON.SETUP'
 319       include 'COMMON.CHAIN'
 320       include 'COMMON.INTERACT'
 321       include 'COMMON.LOCAL'
 322       include 'COMMON.SBRIDGE'
 323       include 'COMMON.TORCNSTR'
 324       include 'COMMON.IOUNITS'
 325       include 'COMMON.DERIV'
 326       include 'COMMON.CONTACTS'
 327       common /przechowalnia/ iturn3_start_all(0:max_fg_procs),
 328      & iturn3_end_all(0:max_fg_procs),iturn4_start_all(0:max_fg_procs),
 329      & iturn4_end_all(0:max_fg_procs),iatel_s_all(0:max_fg_procs),
 330      &iatel_e_all(0:max_fg_procs),ielstart_all(maxres,0:max_fg_procs-1),
 331      & ielend_all(maxres,0:max_fg_procs-1),
 332      & ntask_cont_from_all(0:max_fg_procs-1),
 333      & itask_cont_from_all(0:max_fg_procs-1,0:max_fg_procs-1),
 334      & ntask_cont_to_all(0:max_fg_procs-1),
 335      & itask_cont_to_all(0:max_fg_procs-1,0:max_fg_procs-1)
 336       integer FG_GROUP,CONT_FROM_GROUP,CONT_TO_GROUP
 337       logical scheck,lprint,flag
 338 #ifdef MPI
 339       integer my_sc_int(0:max_fg_Procs-1),my_sc_intt(0:max_fg_Procs),
 340      & my_ele_int(0:max_fg_Procs-1),my_ele_intt(0:max_fg_Procs)
 341 C... Determine the numbers of start and end SC-SC interaction
 342 C... to deal with by current processor.
 343       do i=0,nfgtasks-1
 344         itask_cont_from(i)=fg_rank
 345         itask_cont_to(i)=fg_rank
 346       enddo
 347       lprint=energy_dec
 348       if (lprint)
 349      &write (iout,*) 'INIT_INT_TABLE nres=',nres,' nnt=',nnt,' nct=',nct
 350       n_sc_int_tot=(nct-nnt+1)*(nct-nnt)/2-nss
 351       call int_bounds(n_sc_int_tot,my_sc_inds,my_sc_inde)
 352       if (lprint)
 353      &  write (iout,*) 'Processor',fg_rank,' CG group',kolor,
 354      &  ' absolute rank',MyRank,
 355      &  ' n_sc_int_tot',n_sc_int_tot,' my_sc_inds=',my_sc_inds,
 356      &  ' my_sc_inde',my_sc_inde
 357       ind_sctint=0
 358       iatsc_s=0
 359       iatsc_e=0
 360 #endif
 361 c      lprint=.false.
 362       do i=1,maxres
 363         nint_gr(i)=0
 364         nscp_gr(i)=0
 365         do j=1,maxint_gr
 366           istart(i,1)=0
 367           iend(i,1)=0
 368           ielstart(i)=0
 369           ielend(i)=0
 370           iscpstart(i,1)=0
 371           iscpend(i,1)=0
 372         enddo
 373       enddo
 374       ind_scint=0
 375       ind_scint_old=0
 376 cd    write (iout,*) 'ns=',ns,' nss=',nss,' ihpb,jhpb',
 377 cd   &   (ihpb(i),jhpb(i),i=1,nss)
 378       do i=nnt,nct-1
 379         scheck=.false.
 380         if (dyn_ss) goto 10
 381         do ii=1,nss
 382           if (ihpb(ii).eq.i+nres) then
 383             scheck=.true.
 384             jj=jhpb(ii)-nres
 385             goto 10
 386           endif
 387         enddo
 388    10   continue
 389 cd      write (iout,*) 'i=',i,' scheck=',scheck,' jj=',jj
 390         if (scheck) then
 391           if (jj.eq.i+1) then
 392 #ifdef MPI
 393 c            write (iout,*) 'jj=i+1'
 394             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
 395      & iatsc_s,iatsc_e,i+2,nct,nint_gr(i),istart(i,1),iend(i,1),*12)
 396 #else
 397             nint_gr(i)=1
 398             istart(i,1)=i+2
 399             iend(i,1)=nct
 400 #endif
 401           else if (jj.eq.nct) then
 402 #ifdef MPI
 403 c            write (iout,*) 'jj=nct'
 404             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
 405      &  iatsc_s,iatsc_e,i+1,nct-1,nint_gr(i),istart(i,1),iend(i,1),*12)
 406 #else
 407             nint_gr(i)=1
 408             istart(i,1)=i+1
 409             iend(i,1)=nct-1
 410 #endif
 411           else
 412 #ifdef MPI
 413             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
 414      & iatsc_s,iatsc_e,i+1,jj-1,nint_gr(i),istart(i,1),iend(i,1),*12)
 415             ii=nint_gr(i)+1
 416             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
 417      & iatsc_s,iatsc_e,jj+1,nct,nint_gr(i),istart(i,ii),iend(i,ii),*12)
 418 #else
 419             nint_gr(i)=2
 420             istart(i,1)=i+1
 421             iend(i,1)=jj-1
 422             istart(i,2)=jj+1
 423             iend(i,2)=nct
 424 #endif
 425           endif
 426         else
 427 #ifdef MPI
 428           call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
 429      &    iatsc_s,iatsc_e,i+1,nct,nint_gr(i),istart(i,1),iend(i,1),*12)
 430 #else
 431           nint_gr(i)=1
 432           istart(i,1)=i+1
 433           iend(i,1)=nct
 434           ind_scint=ind_scint+nct-i
 435 #endif
 436         endif
 437 #ifdef MPI
 438         ind_scint_old=ind_scint
 439 #endif
 440       enddo
 441    12 continue
 442 #ifndef MPI
 443       iatsc_s=nnt
 444       iatsc_e=nct-1
 445 #endif
 446       if (iatsc_s.eq.0) iatsc_s=1
 447 #ifdef MPI
 448       if (lprint) write (*,*) 'Processor',fg_rank,' CG Group',kolor,
 449      &   ' absolute rank',myrank,' iatsc_s=',iatsc_s,' iatsc_e=',iatsc_e
 450 #endif
 451       if (lprint) then
 452       write (iout,'(a)') 'Interaction array:'
 453       do i=iatsc_s,iatsc_e
 454         write (iout,'(i3,2(2x,2i3))')
 455      & i,(istart(i,iint),iend(i,iint),iint=1,nint_gr(i))
 456       enddo
 457       endif
 458       ispp=4
 459 #ifdef MPI
 460 C Now partition the electrostatic-interaction array
 461       npept=nct-nnt
 462       nele_int_tot=(npept-ispp)*(npept-ispp+1)/2
 463       call int_bounds(nele_int_tot,my_ele_inds,my_ele_inde)
 464       if (lprint)
 465      & write (*,*) 'Processor',fg_rank,' CG group',kolor,
 466      &  ' absolute rank',MyRank,
 467      &  ' nele_int_tot',nele_int_tot,' my_ele_inds=',my_ele_inds,
 468      &               ' my_ele_inde',my_ele_inde
 469       iatel_s=0
 470       iatel_e=0
 471       ind_eleint=0
 472       ind_eleint_old=0
 473       do i=nnt,nct-3
 474         ijunk=0
 475         call int_partition(ind_eleint,my_ele_inds,my_ele_inde,i,
 476      &    iatel_s,iatel_e,i+ispp,nct-1,ijunk,ielstart(i),ielend(i),*13)
 477       enddo ! i
 478    13 continue
 479       if (iatel_s.eq.0) iatel_s=1
 480       nele_int_tot_vdw=(npept-2)*(npept-2+1)/2
 481 c      write (iout,*) "nele_int_tot_vdw",nele_int_tot_vdw
 482       call int_bounds(nele_int_tot_vdw,my_ele_inds_vdw,my_ele_inde_vdw)
 483 c      write (iout,*) "my_ele_inds_vdw",my_ele_inds_vdw,
 484 c     & " my_ele_inde_vdw",my_ele_inde_vdw
 485       ind_eleint_vdw=0
 486       ind_eleint_vdw_old=0
 487       iatel_s_vdw=0
 488       iatel_e_vdw=0
 489       do i=nnt,nct-3
 490         ijunk=0
 491         call int_partition(ind_eleint_vdw,my_ele_inds_vdw,
 492      &    my_ele_inde_vdw,i,
 493      &    iatel_s_vdw,iatel_e_vdw,i+2,nct-1,ijunk,ielstart_vdw(i),
 494      &    ielend_vdw(i),*15)
 495 c        write (iout,*) i," ielstart_vdw",ielstart_vdw(i),
 496 c     &   " ielend_vdw",ielend_vdw(i)
 497       enddo ! i
 498       if (iatel_s_vdw.eq.0) iatel_s_vdw=1
 499    15 continue
 500 #else
 501       iatel_s=nnt
 502       iatel_e=nct-5
 503       do i=iatel_s,iatel_e
 504         ielstart(i)=i+4
 505         ielend(i)=nct-1
 506       enddo
 507       iatel_s_vdw=nnt
 508       iatel_e_vdw=nct-3
 509       do i=iatel_s_vdw,iatel_e_vdw
 510         ielstart_vdw(i)=i+2
 511         ielend_vdw(i)=nct-1
 512       enddo
 513 #endif
 514       if (lprint) then
 515         write (*,'(a)') 'Processor',fg_rank,' CG group',kolor,
 516      &  ' absolute rank',MyRank
 517         write (iout,*) 'Electrostatic interaction array:'
 518         do i=iatel_s,iatel_e
 519           write (iout,'(i3,2(2x,2i3))') i,ielstart(i),ielend(i)
 520         enddo
 521       endif ! lprint
 522 c     iscp=3
 523       iscp=2
 524 C Partition the SC-p interaction array
 525 #ifdef MPI
 526       nscp_int_tot=(npept-iscp+1)*(npept-iscp+1)
 527       call int_bounds(nscp_int_tot,my_scp_inds,my_scp_inde)
 528       if (lprint) write (iout,*) 'Processor',fg_rank,' CG group',kolor,
 529      &  ' absolute rank',myrank,
 530      &  ' nscp_int_tot',nscp_int_tot,' my_scp_inds=',my_scp_inds,
 531      &               ' my_scp_inde',my_scp_inde
 532       iatscp_s=0
 533       iatscp_e=0
 534       ind_scpint=0
 535       ind_scpint_old=0
 536       do i=nnt,nct-1
 537         if (i.lt.nnt+iscp) then
 538 cd        write (iout,*) 'i.le.nnt+iscp'
 539           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
 540      &      iatscp_s,iatscp_e,i+iscp,nct,nscp_gr(i),iscpstart(i,1),
 541      &      iscpend(i,1),*14)
 542         else if (i.gt.nct-iscp) then
 543 cd        write (iout,*) 'i.gt.nct-iscp'
 544           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
 545      &      iatscp_s,iatscp_e,nnt,i-iscp,nscp_gr(i),iscpstart(i,1),
 546      &      iscpend(i,1),*14)
 547         else
 548           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
 549      &      iatscp_s,iatscp_e,nnt,i-iscp,nscp_gr(i),iscpstart(i,1),
 550      &      iscpend(i,1),*14)
 551           ii=nscp_gr(i)+1
 552           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
 553      &      iatscp_s,iatscp_e,i+iscp,nct,nscp_gr(i),iscpstart(i,ii),
 554      &      iscpend(i,ii),*14)
 555         endif
 556       enddo ! i
 557    14 continue
 558 #else
 559       iatscp_s=nnt
 560       iatscp_e=nct-1
 561       do i=nnt,nct-1
 562         if (i.lt.nnt+iscp) then
 563           nscp_gr(i)=1
 564           iscpstart(i,1)=i+iscp
 565           iscpend(i,1)=nct
 566         elseif (i.gt.nct-iscp) then
 567           nscp_gr(i)=1
 568           iscpstart(i,1)=nnt
 569           iscpend(i,1)=i-iscp
 570         else
 571           nscp_gr(i)=2
 572           iscpstart(i,1)=nnt
 573           iscpend(i,1)=i-iscp
 574           iscpstart(i,2)=i+iscp
 575           iscpend(i,2)=nct
 576         endif
 577       enddo ! i
 578 #endif
 579       if (iatscp_s.eq.0) iatscp_s=1
 580       if (lprint) then
 581         write (iout,'(a)') 'SC-p interaction array:'
 582         do i=iatscp_s,iatscp_e
 583           write (iout,'(i3,2(2x,2i3))')
 584      &         i,(iscpstart(i,j),iscpend(i,j),j=1,nscp_gr(i))
 585         enddo
 586       endif ! lprint
 587 C Partition local interactions
 588 #ifdef MPI
 589       call int_bounds(nres-2,loc_start,loc_end)
 590       loc_start=loc_start+1
 591       loc_end=loc_end+1
 592       call int_bounds(nres-2,ithet_start,ithet_end)
 593       ithet_start=ithet_start+2
 594       ithet_end=ithet_end+2
 595       call int_bounds(nct-nnt-2,iturn3_start,iturn3_end)
 596       iturn3_start=iturn3_start+nnt
 597       iphi_start=iturn3_start+2
 598       iturn3_end=iturn3_end+nnt
 599       iphi_end=iturn3_end+2
 600       iturn3_start=iturn3_start-1
 601       iturn3_end=iturn3_end-1
 602       call int_bounds(nres-3,itau_start,itau_end)
 603       itau_start=itau_start+3
 604       itau_end=itau_end+3
 605       call int_bounds(nres-3,iphi1_start,iphi1_end)
 606       iphi1_start=iphi1_start+3
 607       iphi1_end=iphi1_end+3
 608       call int_bounds(nct-nnt-3,iturn4_start,iturn4_end)
 609       iturn4_start=iturn4_start+nnt
 610       iphid_start=iturn4_start+2
 611       iturn4_end=iturn4_end+nnt
 612       iphid_end=iturn4_end+2
 613       iturn4_start=iturn4_start-1
 614       iturn4_end=iturn4_end-1
 615       call int_bounds(nres-2,ibond_start,ibond_end)
 616       ibond_start=ibond_start+1
 617       ibond_end=ibond_end+1
 618       call int_bounds(nct-nnt,ibondp_start,ibondp_end)
 619       ibondp_start=ibondp_start+nnt
 620       ibondp_end=ibondp_end+nnt
 621       call int_bounds1(nres-1,ivec_start,ivec_end)
 622 c      print *,"Processor",myrank,fg_rank,fg_rank1,
 623 c     &  " ivec_start",ivec_start," ivec_end",ivec_end
 624       iset_start=loc_start+2
 625       iset_end=loc_end+2
 626       if (ndih_constr.eq.0) then
 627         idihconstr_start=1
 628         idihconstr_end=0
 629       else
 630         call int_bounds(ndih_constr,idihconstr_start,idihconstr_end)
 631       endif
 632       if (ntheta_constr.eq.0) then
 633         idihconstr_start=1
 634         idihconstr_end=0
 635       else
 636         call int_bounds
 637      &  (ntheta_constr,ithetaconstr_start,ithetaconstr_end)
 638       endif
 639 c      nsumgrad=(nres-nnt)*(nres-nnt+1)/2
 640 c      nlen=nres-nnt+1
 641       nsumgrad=(nres-nnt)*(nres-nnt+1)/2
 642       nlen=nres-nnt+1
 643       call int_bounds(nsumgrad,ngrad_start,ngrad_end)
 644       igrad_start=((2*nlen+1)
 645      &    -sqrt(float((2*nlen-1)**2-8*(ngrad_start-1))))/2
 646       jgrad_start(igrad_start)=
 647      &    ngrad_start-(2*nlen-igrad_start)*(igrad_start-1)/2
 648      &    +igrad_start
 649       jgrad_end(igrad_start)=nres
 650       igrad_end=((2*nlen+1)
 651      &    -sqrt(float((2*nlen-1)**2-8*(ngrad_end-1))))/2
 652       if (igrad_end.gt.igrad_start) jgrad_start(igrad_end)=igrad_end+1
 653       jgrad_end(igrad_end)=ngrad_end-(2*nlen-igrad_end)*(igrad_end-1)/2
 654      &    +igrad_end
 655       do i=igrad_start+1,igrad_end-1
 656         jgrad_start(i)=i+1
 657         jgrad_end(i)=nres
 658       enddo
 659       if (lprint) then
 660         write (*,*) 'Processor:',fg_rank,' CG group',kolor,
 661      & ' absolute rank',myrank,
 662      & ' loc_start',loc_start,' loc_end',loc_end,
 663      & ' ithet_start',ithet_start,' ithet_end',ithet_end,
 664      & ' iphi_start',iphi_start,' iphi_end',iphi_end,
 665      & ' iphid_start',iphid_start,' iphid_end',iphid_end,
 666      & ' ibond_start',ibond_start,' ibond_end',ibond_end,
 667      & ' ibondp_start',ibondp_start,' ibondp_end',ibondp_end,
 668      & ' iturn3_start',iturn3_start,' iturn3_end',iturn3_end,
 669      & ' iturn4_start',iturn4_start,' iturn4_end',iturn4_end,
 670      & ' ivec_start',ivec_start,' ivec_end',ivec_end,
 671      & ' iset_start',iset_start,' iset_end',iset_end,
 672      & ' idihconstr_start',idihconstr_start,' idihconstr_end',
 673      &   idihconstr_end
 674      & ' ithetaconstr_start',ithetaconstr_start,' ithetaconstr_end',
 675      &   ithetaconstr_end
 676
 677        write (*,*) 'Processor:',fg_rank,myrank,' igrad_start',
 678      &   igrad_start,' igrad_end',igrad_end,' ngrad_start',ngrad_start,
 679      &   ' ngrad_end',ngrad_end
 680        do i=igrad_start,igrad_end
 681          write(*,*) 'Processor:',fg_rank,myrank,i,
 682      &    jgrad_start(i),jgrad_end(i)
 683        enddo
 684       endif
 685       if (nfgtasks.gt.1) then
 686         call MPI_Allgather(ivec_start,1,MPI_INTEGER,ivec_displ(0),1,
 687      &    MPI_INTEGER,FG_COMM1,IERROR)
 688         iaux=ivec_end-ivec_start+1
 689         call MPI_Allgather(iaux,1,MPI_INTEGER,ivec_count(0),1,
 690      &    MPI_INTEGER,FG_COMM1,IERROR)
 691         call MPI_Allgather(iset_start-2,1,MPI_INTEGER,iset_displ(0),1,
 692      &    MPI_INTEGER,FG_COMM,IERROR)
 693         iaux=iset_end-iset_start+1
 694         call MPI_Allgather(iaux,1,MPI_INTEGER,iset_count(0),1,
 695      &    MPI_INTEGER,FG_COMM,IERROR)
 696         call MPI_Allgather(ibond_start,1,MPI_INTEGER,ibond_displ(0),1,
 697      &    MPI_INTEGER,FG_COMM,IERROR)
 698         iaux=ibond_end-ibond_start+1
 699         call MPI_Allgather(iaux,1,MPI_INTEGER,ibond_count(0),1,
 700      &    MPI_INTEGER,FG_COMM,IERROR)
 701         call MPI_Allgather(ithet_start,1,MPI_INTEGER,ithet_displ(0),1,
 702      &    MPI_INTEGER,FG_COMM,IERROR)
 703         iaux=ithet_end-ithet_start+1
 704         call MPI_Allgather(iaux,1,MPI_INTEGER,ithet_count(0),1,
 705      &    MPI_INTEGER,FG_COMM,IERROR)
 706         call MPI_Allgather(iphi_start,1,MPI_INTEGER,iphi_displ(0),1,
 707      &    MPI_INTEGER,FG_COMM,IERROR)
 708         iaux=iphi_end-iphi_start+1
 709         call MPI_Allgather(iaux,1,MPI_INTEGER,iphi_count(0),1,
 710      &    MPI_INTEGER,FG_COMM,IERROR)
 711         call MPI_Allgather(iphi1_start,1,MPI_INTEGER,iphi1_displ(0),1,
 712      &    MPI_INTEGER,FG_COMM,IERROR)
 713         iaux=iphi1_end-iphi1_start+1
 714         call MPI_Allgather(iaux,1,MPI_INTEGER,iphi1_count(0),1,
 715      &    MPI_INTEGER,FG_COMM,IERROR)
 716         do i=0,maxprocs-1
 717           do j=1,maxres
 718             ielstart_all(j,i)=0
 719             ielend_all(j,i)=0
 720           enddo
 721         enddo
 722         call MPI_Allgather(iturn3_start,1,MPI_INTEGER,
 723      &    iturn3_start_all(0),1,MPI_INTEGER,FG_COMM,IERROR)
 724         call MPI_Allgather(iturn4_start,1,MPI_INTEGER,
 725      &    iturn4_start_all(0),1,MPI_INTEGER,FG_COMM,IERROR)
 726         call MPI_Allgather(iturn3_end,1,MPI_INTEGER,
 727      &    iturn3_end_all(0),1,MPI_INTEGER,FG_COMM,IERROR)
 728         call MPI_Allgather(iturn4_end,1,MPI_INTEGER,
 729      &    iturn4_end_all(0),1,MPI_INTEGER,FG_COMM,IERROR)
 730         call MPI_Allgather(iatel_s,1,MPI_INTEGER,
 731      &    iatel_s_all(0),1,MPI_INTEGER,FG_COMM,IERROR)
 732         call MPI_Allgather(iatel_e,1,MPI_INTEGER,
 733      &    iatel_e_all(0),1,MPI_INTEGER,FG_COMM,IERROR)
 734         call MPI_Allgather(ielstart(1),maxres,MPI_INTEGER,
 735      &    ielstart_all(1,0),maxres,MPI_INTEGER,FG_COMM,IERROR)
 736         call MPI_Allgather(ielend(1),maxres,MPI_INTEGER,
 737      &    ielend_all(1,0),maxres,MPI_INTEGER,FG_COMM,IERROR)
 738         if (lprint) then
 739         write (iout,*) "iatel_s_all",(iatel_s_all(i),i=0,nfgtasks)
 740         write (iout,*) "iatel_e_all",(iatel_e_all(i),i=0,nfgtasks)
 741         write (iout,*) "iturn3_start_all",
 742      &    (iturn3_start_all(i),i=0,nfgtasks-1)
 743         write (iout,*) "iturn3_end_all",
 744      &    (iturn3_end_all(i),i=0,nfgtasks-1)
 745         write (iout,*) "iturn4_start_all",
 746      &    (iturn4_start_all(i),i=0,nfgtasks-1)
 747         write (iout,*) "iturn4_end_all",
 748      &    (iturn4_end_all(i),i=0,nfgtasks-1)
 749         write (iout,*) "The ielstart_all array"
 750         do i=nnt,nct
 751           write (iout,'(20i4)') i,(ielstart_all(i,j),j=0,nfgtasks-1)
 752         enddo
 753         write (iout,*) "The ielend_all array"
 754         do i=nnt,nct
 755           write (iout,'(20i4)') i,(ielend_all(i,j),j=0,nfgtasks-1)
 756         enddo
 757         call flush(iout)
 758         endif
 759         ntask_cont_from=0
 760         ntask_cont_to=0
 761         itask_cont_from(0)=fg_rank
 762         itask_cont_to(0)=fg_rank
 763         flag=.false.
 764         do ii=iturn3_start,iturn3_end
 765           call add_int(ii,ii+2,iturn3_sent(1,ii),
 766      &                 ntask_cont_to,itask_cont_to,flag)
 767         enddo
 768         do ii=iturn4_start,iturn4_end
 769           call add_int(ii,ii+3,iturn4_sent(1,ii),
 770      &                 ntask_cont_to,itask_cont_to,flag)
 771         enddo
 772         do ii=iturn3_start,iturn3_end
 773           call add_int_from(ii,ii+2,ntask_cont_from,itask_cont_from)
 774         enddo
 775         do ii=iturn4_start,iturn4_end
 776           call add_int_from(ii,ii+3,ntask_cont_from,itask_cont_from)
 777         enddo
 778         if (lprint) then
 779         write (iout,*) "After turn3 ntask_cont_from",ntask_cont_from,
 780      &   " ntask_cont_to",ntask_cont_to
 781         write (iout,*) "itask_cont_from",
 782      &    (itask_cont_from(i),i=1,ntask_cont_from)
 783         write (iout,*) "itask_cont_to",
 784      &    (itask_cont_to(i),i=1,ntask_cont_to)
 785         call flush(iout)
 786         endif
 787 c        write (iout,*) "Loop forward"
 788 c        call flush(iout)
 789         do i=iatel_s,iatel_e
 790 c          write (iout,*) "from loop i=",i
 791 c          call flush(iout)
 792           do j=ielstart(i),ielend(i)
 793             call add_int_from(i,j,ntask_cont_from,itask_cont_from)
 794           enddo
 795         enddo
 796 c        write (iout,*) "Loop backward iatel_e-1",iatel_e-1,
 797 c     &     " iatel_e",iatel_e
 798 c        call flush(iout)
 799         nat_sent=0
 800         do i=iatel_s,iatel_e
 801 c          write (iout,*) "i",i," ielstart",ielstart(i),
 802 c     &      " ielend",ielend(i)
 803 c          call flush(iout)
 804           flag=.false.
 805           do j=ielstart(i),ielend(i)
 806             call add_int(i,j,iint_sent(1,j,nat_sent+1),ntask_cont_to,
 807      &                   itask_cont_to,flag)
 808           enddo
 809           if (flag) then
 810             nat_sent=nat_sent+1
 811             iat_sent(nat_sent)=i
 812           endif
 813         enddo
 814         if (lprint) then
 815         write (iout,*)"After longrange ntask_cont_from",ntask_cont_from,
 816      &   " ntask_cont_to",ntask_cont_to
 817         write (iout,*) "itask_cont_from",
 818      &    (itask_cont_from(i),i=1,ntask_cont_from)
 819         write (iout,*) "itask_cont_to",
 820      &    (itask_cont_to(i),i=1,ntask_cont_to)
 821         call flush(iout)
 822         write (iout,*) "iint_sent"
 823         do i=1,nat_sent
 824           ii=iat_sent(i)
 825           write (iout,'(20i4)') ii,(j,(iint_sent(k,j,i),k=1,4),
 826      &      j=ielstart(ii),ielend(ii))
 827         enddo
 828         write (iout,*) "iturn3_sent iturn3_start",iturn3_start,
 829      &    " iturn3_end",iturn3_end
 830         write (iout,'(20i4)') (i,(iturn3_sent(j,i),j=1,4),
 831      &      i=iturn3_start,iturn3_end)
 832         write (iout,*) "iturn4_sent iturn4_start",iturn4_start,
 833      &    " iturn4_end",iturn4_end
 834         write (iout,'(20i4)') (i,(iturn4_sent(j,i),j=1,4),
 835      &      i=iturn4_start,iturn4_end)
 836         call flush(iout)
 837         endif
 838         call MPI_Gather(ntask_cont_from,1,MPI_INTEGER,
 839      &   ntask_cont_from_all,1,MPI_INTEGER,king,FG_COMM,IERR)
 840 c        write (iout,*) "Gather ntask_cont_from ended"
 841 c        call flush(iout)
 842         call MPI_Gather(itask_cont_from(0),max_fg_procs,MPI_INTEGER,
 843      &   itask_cont_from_all(0,0),max_fg_procs,MPI_INTEGER,king,
 844      &   FG_COMM,IERR)
 845 c        write (iout,*) "Gather itask_cont_from ended"
 846 c        call flush(iout)
 847         call MPI_Gather(ntask_cont_to,1,MPI_INTEGER,ntask_cont_to_all,
 848      &   1,MPI_INTEGER,king,FG_COMM,IERR)
 849 c        write (iout,*) "Gather ntask_cont_to ended"
 850 c        call flush(iout)
 851         call MPI_Gather(itask_cont_to,max_fg_procs,MPI_INTEGER,
 852      &   itask_cont_to_all,max_fg_procs,MPI_INTEGER,king,FG_COMM,IERR)
 853 c        write (iout,*) "Gather itask_cont_to ended"
 854 c        call flush(iout)
 855         if (fg_rank.eq.king) then
 856           write (iout,*)"Contact receive task map (proc, #tasks, tasks)"
 857           do i=0,nfgtasks-1
 858             write (iout,'(20i4)') i,ntask_cont_from_all(i),
 859      &       (itask_cont_from_all(j,i),j=1,ntask_cont_from_all(i))
 860           enddo
 861           write (iout,*)
 862           call flush(iout)
 863           write (iout,*) "Contact send task map (proc, #tasks, tasks)"
 864           do i=0,nfgtasks-1
 865             write (iout,'(20i4)') i,ntask_cont_to_all(i),
 866      &       (itask_cont_to_all(j,i),j=1,ntask_cont_to_all(i))
 867           enddo
 868           write (iout,*)
 869           call flush(iout)
 870 C Check if every send will have a matching receive
 871           ncheck_to=0
 872           ncheck_from=0
 873           do i=0,nfgtasks-1
 874             ncheck_to=ncheck_to+ntask_cont_to_all(i)
 875             ncheck_from=ncheck_from+ntask_cont_from_all(i)
 876           enddo
 877           write (iout,*) "Control sums",ncheck_from,ncheck_to
 878           if (ncheck_from.ne.ncheck_to) then
 879             write (iout,*) "Error: #receive differs from #send."
 880             write (iout,*) "Terminating program...!"
 881             call flush(iout)
 882             flag=.false.
 883           else
 884             flag=.true.
 885             do i=0,nfgtasks-1
 886               do j=1,ntask_cont_to_all(i)
 887                 ii=itask_cont_to_all(j,i)
 888                 do k=1,ntask_cont_from_all(ii)
 889                   if (itask_cont_from_all(k,ii).eq.i) then
 890                     if(lprint)write(iout,*)"Matching send/receive",i,ii
 891                     exit
 892                   endif
 893                 enddo
 894                 if (k.eq.ntask_cont_from_all(ii)+1) then
 895                   flag=.false.
 896                   write (iout,*) "Error: send by",j," to",ii,
 897      &            " would have no matching receive"
 898                 endif
 899               enddo
 900             enddo
 901           endif
 902           if (.not.flag) then
 903             write (iout,*) "Unmatched sends; terminating program"
 904             call flush(iout)
 905           endif
 906         endif
 907         call MPI_Bcast(flag,1,MPI_LOGICAL,king,FG_COMM,IERROR)
 908 c        write (iout,*) "flag broadcast ended flag=",flag
 909 c        call flush(iout)
 910         if (.not.flag) then
 911           call MPI_Finalize(IERROR)
 912           stop "Error in INIT_INT_TABLE: unmatched send/receive."
 913         endif
 914         call MPI_Comm_group(FG_COMM,fg_group,IERR)
 915 c        write (iout,*) "MPI_Comm_group ended"
 916 c        call flush(iout)
 917         call MPI_Group_incl(fg_group,ntask_cont_from+1,
 918      &    itask_cont_from(0),CONT_FROM_GROUP,IERR)
 919         call MPI_Group_incl(fg_group,ntask_cont_to+1,itask_cont_to(0),
 920      &    CONT_TO_GROUP,IERR)
 921         do i=1,nat_sent
 922           ii=iat_sent(i)
 923           iaux=4*(ielend(ii)-ielstart(ii)+1)
 924           call MPI_Group_translate_ranks(fg_group,iaux,
 925      &      iint_sent(1,ielstart(ii),i),CONT_TO_GROUP,
 926      &      iint_sent_local(1,ielstart(ii),i),IERR )
 927 c          write (iout,*) "Ranks translated i=",i
 928 c          call flush(iout)
 929         enddo
 930         iaux=4*(iturn3_end-iturn3_start+1)
 931         call MPI_Group_translate_ranks(fg_group,iaux,
 932      &      iturn3_sent(1,iturn3_start),CONT_TO_GROUP,
 933      &      iturn3_sent_local(1,iturn3_start),IERR)
 934         iaux=4*(iturn4_end-iturn4_start+1)
 935         call MPI_Group_translate_ranks(fg_group,iaux,
 936      &      iturn4_sent(1,iturn4_start),CONT_TO_GROUP,
 937      &      iturn4_sent_local(1,iturn4_start),IERR)
 938         if (lprint) then
 939         write (iout,*) "iint_sent_local"
 940         do i=1,nat_sent
 941           ii=iat_sent(i)
 942           write (iout,'(20i4)') ii,(j,(iint_sent_local(k,j,i),k=1,4),
 943      &      j=ielstart(ii),ielend(ii))
 944           call flush(iout)
 945         enddo
 946         write (iout,*) "iturn3_sent_local iturn3_start",iturn3_start,
 947      &    " iturn3_end",iturn3_end
 948         write (iout,'(20i4)') (i,(iturn3_sent_local(j,i),j=1,4),
 949      &      i=iturn3_start,iturn3_end)
 950         write (iout,*) "iturn4_sent_local iturn4_start",iturn4_start,
 951      &    " iturn4_end",iturn4_end
 952         write (iout,'(20i4)') (i,(iturn4_sent_local(j,i),j=1,4),
 953      &      i=iturn4_start,iturn4_end)
 954         call flush(iout)
 955         endif
 956         call MPI_Group_free(fg_group,ierr)
 957         call MPI_Group_free(cont_from_group,ierr)
 958         call MPI_Group_free(cont_to_group,ierr)
 959         call MPI_Type_contiguous(3,MPI_DOUBLE_PRECISION,MPI_UYZ,IERROR)
 960         call MPI_Type_commit(MPI_UYZ,IERROR)
 961         call MPI_Type_contiguous(18,MPI_DOUBLE_PRECISION,MPI_UYZGRAD,
 962      &    IERROR)
 963         call MPI_Type_commit(MPI_UYZGRAD,IERROR)
 964         call MPI_Type_contiguous(2,MPI_DOUBLE_PRECISION,MPI_MU,IERROR)
 965         call MPI_Type_commit(MPI_MU,IERROR)
 966         call MPI_Type_contiguous(4,MPI_DOUBLE_PRECISION,MPI_MAT1,IERROR)
 967         call MPI_Type_commit(MPI_MAT1,IERROR)
 968         call MPI_Type_contiguous(8,MPI_DOUBLE_PRECISION,MPI_MAT2,IERROR)
 969         call MPI_Type_commit(MPI_MAT2,IERROR)
 970         call MPI_Type_contiguous(6,MPI_DOUBLE_PRECISION,MPI_THET,IERROR)
 971         call MPI_Type_commit(MPI_THET,IERROR)
 972         call MPI_Type_contiguous(9,MPI_DOUBLE_PRECISION,MPI_GAM,IERROR)
 973         call MPI_Type_commit(MPI_GAM,IERROR)
 974 #ifndef MATGATHER
 975 c 9/22/08 Derived types to send matrices which appear in correlation terms
 976         do i=0,nfgtasks-1
 977           if (ivec_count(i).eq.ivec_count(0)) then
 978             lentyp(i)=0
 979           else
 980             lentyp(i)=1
 981           endif
 982         enddo
 983         do ind_typ=lentyp(0),lentyp(nfgtasks-1)
 984         if (ind_typ.eq.0) then
 985           ichunk=ivec_count(0)
 986         else
 987           ichunk=ivec_count(1)
 988         endif
 989 c        do i=1,4
 990 c          blocklengths(i)=4
 991 c        enddo
 992 c        displs(1)=0
 993 c        do i=2,4
 994 c          displs(i)=displs(i-1)+blocklengths(i-1)*maxres
 995 c        enddo
 996 c        do i=1,4
 997 c          blocklengths(i)=blocklengths(i)*ichunk
 998 c        enddo
 999 c        write (iout,*) "blocklengths and displs"
1000 c        do i=1,4
1001 c          write (iout,*) i,blocklengths(i),displs(i)
1002 c        enddo
1003 c        call flush(iout)
1004 c        call MPI_Type_indexed(4,blocklengths(1),displs(1),
1005 c     &    MPI_DOUBLE_PRECISION,MPI_ROTAT1(ind_typ),IERROR)
1006 c        call MPI_Type_commit(MPI_ROTAT1(ind_typ),IERROR)
1007 c        write (iout,*) "MPI_ROTAT1",MPI_ROTAT1
1008 c        do i=1,4
1009 c          blocklengths(i)=2
1010 c        enddo
1011 c        displs(1)=0
1012 c        do i=2,4
1013 c          displs(i)=displs(i-1)+blocklengths(i-1)*maxres
1014 c        enddo
1015 c        do i=1,4
1016 c          blocklengths(i)=blocklengths(i)*ichunk
1017 c        enddo
1018 c        write (iout,*) "blocklengths and displs"
1019 c        do i=1,4
1020 c          write (iout,*) i,blocklengths(i),displs(i)
1021 c        enddo
1022 c        call flush(iout)
1023 c        call MPI_Type_indexed(4,blocklengths(1),displs(1),
1024 c     &    MPI_DOUBLE_PRECISION,MPI_ROTAT2(ind_typ),IERROR)
1025 c        call MPI_Type_commit(MPI_ROTAT2(ind_typ),IERROR)
1026 c        write (iout,*) "MPI_ROTAT2",MPI_ROTAT2
1027         do i=1,8
1028           blocklengths(i)=2
1029         enddo
1030         displs(1)=0
1031         do i=2,8
1032           displs(i)=displs(i-1)+blocklengths(i-1)*maxres
1033         enddo
1034         do i=1,15
1035           blocklengths(i)=blocklengths(i)*ichunk
1036         enddo
1037         call MPI_Type_indexed(8,blocklengths,displs,
1038      &    MPI_DOUBLE_PRECISION,MPI_PRECOMP11(ind_typ),IERROR)
1039         call MPI_Type_commit(MPI_PRECOMP11(ind_typ),IERROR)
1040         do i=1,8
1041           blocklengths(i)=4
1042         enddo
1043         displs(1)=0
1044         do i=2,8
1045           displs(i)=displs(i-1)+blocklengths(i-1)*maxres
1046         enddo
1047         do i=1,15
1048           blocklengths(i)=blocklengths(i)*ichunk
1049         enddo
1050         call MPI_Type_indexed(8,blocklengths,displs,
1051      &    MPI_DOUBLE_PRECISION,MPI_PRECOMP12(ind_typ),IERROR)
1052         call MPI_Type_commit(MPI_PRECOMP12(ind_typ),IERROR)
1053         do i=1,6
1054           blocklengths(i)=4
1055         enddo
1056         displs(1)=0
1057         do i=2,6
1058           displs(i)=displs(i-1)+blocklengths(i-1)*maxres
1059         enddo
1060         do i=1,6
1061           blocklengths(i)=blocklengths(i)*ichunk
1062         enddo
1063         call MPI_Type_indexed(6,blocklengths,displs,
1064      &    MPI_DOUBLE_PRECISION,MPI_PRECOMP22(ind_typ),IERROR)
1065         call MPI_Type_commit(MPI_PRECOMP22(ind_typ),IERROR)
1066         do i=1,2
1067           blocklengths(i)=8
1068         enddo
1069         displs(1)=0
1070         do i=2,2
1071           displs(i)=displs(i-1)+blocklengths(i-1)*maxres
1072         enddo
1073         do i=1,2
1074           blocklengths(i)=blocklengths(i)*ichunk
1075         enddo
1076         call MPI_Type_indexed(2,blocklengths,displs,
1077      &    MPI_DOUBLE_PRECISION,MPI_PRECOMP23(ind_typ),IERROR)
1078         call MPI_Type_commit(MPI_PRECOMP23(ind_typ),IERROR)
1079         do i=1,4
1080           blocklengths(i)=1
1081         enddo
1082         displs(1)=0
1083         do i=2,4
1084           displs(i)=displs(i-1)+blocklengths(i-1)*maxres
1085         enddo
1086         do i=1,4
1087           blocklengths(i)=blocklengths(i)*ichunk
1088         enddo
1089         call MPI_Type_indexed(4,blocklengths,displs,
1090      &    MPI_DOUBLE_PRECISION,MPI_ROTAT_OLD(ind_typ),IERROR)
1091         call MPI_Type_commit(MPI_ROTAT_OLD(ind_typ),IERROR)
1092         enddo
1093 #endif
1094       endif
1095       iint_start=ivec_start+1
1096       iint_end=ivec_end+1
1097       do i=0,nfgtasks-1
1098           iint_count(i)=ivec_count(i)
1099           iint_displ(i)=ivec_displ(i)
1100           ivec_displ(i)=ivec_displ(i)-1
1101           iset_displ(i)=iset_displ(i)-1
1102           ithet_displ(i)=ithet_displ(i)-1
1103           iphi_displ(i)=iphi_displ(i)-1
1104           iphi1_displ(i)=iphi1_displ(i)-1
1105           ibond_displ(i)=ibond_displ(i)-1
1106       enddo
1107       if (nfgtasks.gt.1 .and. fg_rank.eq.king
1108      &     .and. (me.eq.0 .or. .not. out1file)) then
1109         write (iout,*) "IVEC_DISPL, IVEC_COUNT, ISET_START, ISET_COUNT"
1110         do i=0,nfgtasks-1
1111           write (iout,*) i,ivec_displ(i),ivec_count(i),iset_displ(i),
1112      &      iset_count(i)
1113         enddo
1114         write (iout,*) "iphi_start",iphi_start," iphi_end",iphi_end,
1115      &    " iphi1_start",iphi1_start," iphi1_end",iphi1_end
1116         write (iout,*)"IPHI_COUNT, IPHI_DISPL, IPHI1_COUNT, IPHI1_DISPL"
1117         do i=0,nfgtasks-1
1118           write (iout,*) i,iphi_count(i),iphi_displ(i),iphi1_count(i),
1119      &      iphi1_displ(i)
1120         enddo
1121         write(iout,'(i10,a,i10,a,i10,a/a,i3,a)') n_sc_int_tot,' SC-SC ',
1122      & nele_int_tot,' electrostatic and ',nscp_int_tot,
1123      & ' SC-p interactions','were distributed among',nfgtasks,
1124      & ' fine-grain processors.'
1125       endif
1126 #else
1127       loc_start=2
1128       loc_end=nres-1
1129       ithet_start=3
1130       ithet_end=nres
1131       iturn3_start=nnt
1132       iturn3_end=nct-3
1133       iturn4_start=nnt
1134       iturn4_end=nct-4
1135       iphi_start=nnt+3
1136       iphi_end=nct
1137       iphi1_start=4
1138       iphi1_end=nres
1139       idihconstr_start=1
1140       idihconstr_end=ndih_constr
1141       ithetaconstr_start=1
1142       ithetaconstr_end=ntheta_constr
1143       iphid_start=iphi_start
1144       iphid_end=iphi_end-1
1145       itau_start=4
1146       itau_end=nres
1147       ibond_start=2
1148       ibond_end=nres-1
1149       ibondp_start=nnt
1150       ibondp_end=nct-1
1151       ivec_start=1
1152       ivec_end=nres-1
1153       iset_start=3
1154       iset_end=nres+1
1155       iint_start=2
1156       iint_end=nres-1
1157 #endif
1158       return
1159       end
1160 #ifdef MPI
1161 c---------------------------------------------------------------------------
1162       subroutine add_int(ii,jj,itask,ntask_cont_to,itask_cont_to,flag)
1163       implicit none
1164       include "DIMENSIONS"
1165       include "COMMON.INTERACT"
1166       include "COMMON.SETUP"
1167       include "COMMON.IOUNITS"
1168       integer ii,jj,itask(4),ntask_cont_to,
1169      &itask_cont_to(0:max_fg_procs-1)
1170       logical flag
1171       integer iturn3_start_all,iturn3_end_all,iturn4_start_all,
1172      & iturn4_end_all,iatel_s_all,iatel_e_all,ielstart_all,ielend_all
1173       common /przechowalnia/ iturn3_start_all(0:max_fg_procs),
1174      & iturn3_end_all(0:max_fg_procs),iturn4_start_all(0:max_fg_procs),
1175      & iturn4_end_all(0:max_fg_procs),iatel_s_all(0:max_fg_procs),
1176      &iatel_e_all(0:max_fg_procs),ielstart_all(maxres,0:max_fg_procs-1),
1177      & ielend_all(maxres,0:max_fg_procs-1)
1178       integer iproc,isent,k,l
1179 c Determines whether to send interaction ii,jj to other processors; a given
1180 c interaction can be sent to at most 2 processors.
1181 c Sets flag=.true. if interaction ii,jj needs to be sent to at least
1182 c one processor, otherwise flag is unchanged from the input value.
1183       isent=0
1184       itask(1)=fg_rank
1185       itask(2)=fg_rank
1186       itask(3)=fg_rank
1187       itask(4)=fg_rank
1188 c      write (iout,*) "ii",ii," jj",jj
1189 c Loop over processors to check if anybody could need interaction ii,jj
1190       do iproc=0,fg_rank-1
1191 c Check if the interaction matches any turn3 at iproc
1192         do k=iturn3_start_all(iproc),iturn3_end_all(iproc)
1193           l=k+2
1194           if (k.eq.ii-1 .and. l.eq.jj-1 .or. k.eq.ii-1 .and. l.eq.jj+1
1195      &   .or. k.eq.ii+1 .and. l.eq.jj+1 .or. k.eq.ii+1 .and. l.eq.jj-1)
1196      &    then
1197 c            write (iout,*) "turn3 to iproc",iproc," ij",ii,jj,"kl",k,l
1198 c            call flush(iout)
1199             flag=.true.
1200             if (iproc.ne.itask(1).and.iproc.ne.itask(2)
1201      &        .and.iproc.ne.itask(3).and.iproc.ne.itask(4)) then
1202               isent=isent+1
1203               itask(isent)=iproc
1204               call add_task(iproc,ntask_cont_to,itask_cont_to)
1205             endif
1206           endif
1207         enddo
1208 C Check if the interaction matches any turn4 at iproc
1209         do k=iturn4_start_all(iproc),iturn4_end_all(iproc)
1210           l=k+3
1211           if (k.eq.ii-1 .and. l.eq.jj-1 .or. k.eq.ii-1 .and. l.eq.jj+1
1212      &   .or. k.eq.ii+1 .and. l.eq.jj+1 .or. k.eq.ii+1 .and. l.eq.jj-1)
1213      &    then
1214 c            write (iout,*) "turn3 to iproc",iproc," ij",ii,jj," kl",k,l
1215 c            call flush(iout)
1216             flag=.true.
1217             if (iproc.ne.itask(1).and.iproc.ne.itask(2)
1218      &        .and.iproc.ne.itask(3).and.iproc.ne.itask(4)) then
1219               isent=isent+1
1220               itask(isent)=iproc
1221               call add_task(iproc,ntask_cont_to,itask_cont_to)
1222             endif
1223           endif
1224         enddo
1225         if (iatel_s_all(iproc).gt.0 .and. iatel_e_all(iproc).gt.0 .and.
1226      &  iatel_s_all(iproc).le.ii-1 .and. iatel_e_all(iproc).ge.ii-1)then
1227           if (ielstart_all(ii-1,iproc).le.jj-1.and.
1228      &        ielend_all(ii-1,iproc).ge.jj-1) then
1229             flag=.true.
1230             if (iproc.ne.itask(1).and.iproc.ne.itask(2)
1231      &        .and.iproc.ne.itask(3).and.iproc.ne.itask(4)) then
1232               isent=isent+1
1233               itask(isent)=iproc
1234               call add_task(iproc,ntask_cont_to,itask_cont_to)
1235             endif
1236           endif
1237           if (ielstart_all(ii-1,iproc).le.jj+1.and.
1238      &        ielend_all(ii-1,iproc).ge.jj+1) then
1239             flag=.true.
1240             if (iproc.ne.itask(1).and.iproc.ne.itask(2)
1241      &        .and.iproc.ne.itask(3).and.iproc.ne.itask(4)) then
1242               isent=isent+1
1243               itask(isent)=iproc
1244               call add_task(iproc,ntask_cont_to,itask_cont_to)
1245             endif
1246           endif
1247         endif
1248       enddo
1249       return
1250       end
1251 c---------------------------------------------------------------------------
1252       subroutine add_int_from(ii,jj,ntask_cont_from,itask_cont_from)
1253       implicit none
1254       include "DIMENSIONS"
1255       include "COMMON.INTERACT"
1256       include "COMMON.SETUP"
1257       include "COMMON.IOUNITS"
1258       integer ii,jj,itask(2),ntask_cont_from,
1259      & itask_cont_from(0:max_fg_procs-1)
1260       logical flag
1261       integer iturn3_start_all,iturn3_end_all,iturn4_start_all,
1262      & iturn4_end_all,iatel_s_all,iatel_e_all,ielstart_all,ielend_all
1263       common /przechowalnia/ iturn3_start_all(0:max_fg_procs),
1264      & iturn3_end_all(0:max_fg_procs),iturn4_start_all(0:max_fg_procs),
1265      & iturn4_end_all(0:max_fg_procs),iatel_s_all(0:max_fg_procs),
1266      &iatel_e_all(0:max_fg_procs),ielstart_all(maxres,0:max_fg_procs-1),
1267      & ielend_all(maxres,0:max_fg_procs-1)
1268       integer iproc,k,l
1269       do iproc=fg_rank+1,nfgtasks-1
1270         do k=iturn3_start_all(iproc),iturn3_end_all(iproc)
1271           l=k+2
1272           if (k.eq.ii+1 .and. l.eq.jj+1 .or. k.eq.ii+1.and.l.eq.jj-1
1273      &   .or. k.eq.ii-1 .and. l.eq.jj-1 .or. k.eq.ii-1 .and. l.eq.jj+1)
1274      &    then
1275 c            write (iout,*)"turn3 from iproc",iproc," ij",ii,jj," kl",k,l
1276             call add_task(iproc,ntask_cont_from,itask_cont_from)
1277           endif
1278         enddo
1279         do k=iturn4_start_all(iproc),iturn4_end_all(iproc)
1280           l=k+3
1281           if (k.eq.ii+1 .and. l.eq.jj+1 .or. k.eq.ii+1.and.l.eq.jj-1
1282      &   .or. k.eq.ii-1 .and. l.eq.jj-1 .or. k.eq.ii-1 .and. l.eq.jj+1)
1283      &    then
1284 c            write (iout,*)"turn4 from iproc",iproc," ij",ii,jj," kl",k,l
1285             call add_task(iproc,ntask_cont_from,itask_cont_from)
1286           endif
1287         enddo
1288         if (iatel_s_all(iproc).gt.0 .and. iatel_e_all(iproc).gt.0) then
1289           if (ii+1.ge.iatel_s_all(iproc).and.ii+1.le.iatel_e_all(iproc))
1290      &    then
1291             if (jj+1.ge.ielstart_all(ii+1,iproc).and.
1292      &          jj+1.le.ielend_all(ii+1,iproc)) then
1293               call add_task(iproc,ntask_cont_from,itask_cont_from)
1294             endif
1295             if (jj-1.ge.ielstart_all(ii+1,iproc).and.
1296      &          jj-1.le.ielend_all(ii+1,iproc)) then
1297               call add_task(iproc,ntask_cont_from,itask_cont_from)
1298             endif
1299           endif
1300           if (ii-1.ge.iatel_s_all(iproc).and.ii-1.le.iatel_e_all(iproc))
1301      &    then
1302             if (jj-1.ge.ielstart_all(ii-1,iproc).and.
1303      &          jj-1.le.ielend_all(ii-1,iproc)) then
1304               call add_task(iproc,ntask_cont_from,itask_cont_from)
1305             endif
1306             if (jj+1.ge.ielstart_all(ii-1,iproc).and.
1307      &          jj+1.le.ielend_all(ii-1,iproc)) then
1308                call add_task(iproc,ntask_cont_from,itask_cont_from)
1309             endif
1310           endif
1311         endif
1312       enddo
1313       return
1314       end
1315 c---------------------------------------------------------------------------
1316       subroutine add_task(iproc,ntask_cont,itask_cont)
1317       implicit none
1318       include "DIMENSIONS"
1319       integer iproc,ntask_cont,itask_cont(0:max_fg_procs-1)
1320       integer ii
1321       do ii=1,ntask_cont
1322         if (itask_cont(ii).eq.iproc) return
1323       enddo
1324       ntask_cont=ntask_cont+1
1325       itask_cont(ntask_cont)=iproc
1326       return
1327       end
1328 c---------------------------------------------------------------------------
1329       subroutine int_bounds(total_ints,lower_bound,upper_bound)
1330       implicit real*8 (a-h,o-z)
1331       include 'DIMENSIONS'
1332       include 'mpif.h'
1333       include 'COMMON.SETUP'
1334       integer total_ints,lower_bound,upper_bound
1335       integer int4proc(0:max_fg_procs),sint4proc(0:max_fg_procs)
1336       nint=total_ints/nfgtasks
1337       do i=1,nfgtasks
1338         int4proc(i-1)=nint
1339       enddo
1340       nexcess=total_ints-nint*nfgtasks
1341       do i=1,nexcess
1342         int4proc(nfgtasks-i)=int4proc(nfgtasks-i)+1
1343       enddo
1344       lower_bound=0
1345       do i=0,fg_rank-1
1346         lower_bound=lower_bound+int4proc(i)
1347       enddo
1348       upper_bound=lower_bound+int4proc(fg_rank)
1349       lower_bound=lower_bound+1
1350       return
1351       end
1352 c---------------------------------------------------------------------------
1353       subroutine int_bounds1(total_ints,lower_bound,upper_bound)
1354       implicit real*8 (a-h,o-z)
1355       include 'DIMENSIONS'
1356       include 'mpif.h'
1357       include 'COMMON.SETUP'
1358       integer total_ints,lower_bound,upper_bound
1359       integer int4proc(0:max_fg_procs),sint4proc(0:max_fg_procs)
1360       nint=total_ints/nfgtasks1
1361       do i=1,nfgtasks1
1362         int4proc(i-1)=nint
1363       enddo
1364       nexcess=total_ints-nint*nfgtasks1
1365       do i=1,nexcess
1366         int4proc(nfgtasks1-i)=int4proc(nfgtasks1-i)+1
1367       enddo
1368       lower_bound=0
1369       do i=0,fg_rank1-1
1370         lower_bound=lower_bound+int4proc(i)
1371       enddo
1372       upper_bound=lower_bound+int4proc(fg_rank1)
1373       lower_bound=lower_bound+1
1374       return
1375       end
1376 c---------------------------------------------------------------------------
1377       subroutine int_partition(int_index,lower_index,upper_index,atom,
1378      & at_start,at_end,first_atom,last_atom,int_gr,jat_start,jat_end,*)
1379       implicit real*8 (a-h,o-z)
1380       include 'DIMENSIONS'
1381       include 'COMMON.IOUNITS'
1382       integer int_index,lower_index,upper_index,atom,at_start,at_end,
1383      & first_atom,last_atom,int_gr,jat_start,jat_end
1384       logical lprn
1385       lprn=.false.
1386       if (lprn) write (iout,*) 'int_index=',int_index
1387       int_index_old=int_index
1388       int_index=int_index+last_atom-first_atom+1
1389       if (lprn)
1390      &   write (iout,*) 'int_index=',int_index,
1391      &               ' int_index_old',int_index_old,
1392      &               ' lower_index=',lower_index,
1393      &               ' upper_index=',upper_index,
1394      &               ' atom=',atom,' first_atom=',first_atom,
1395      &               ' last_atom=',last_atom
1396       if (int_index.ge.lower_index) then
1397         int_gr=int_gr+1
1398         if (at_start.eq.0) then
1399           at_start=atom
1400           jat_start=first_atom-1+lower_index-int_index_old
1401         else
1402           jat_start=first_atom
1403         endif
1404         if (lprn) write (iout,*) 'jat_start',jat_start
1405         if (int_index.ge.upper_index) then
1406           at_end=atom
1407           jat_end=first_atom-1+upper_index-int_index_old
1408           return1
1409         else
1410           jat_end=last_atom
1411         endif
1412         if (lprn) write (iout,*) 'jat_end',jat_end
1413       endif
1414       return
1415       end
1416 #endif
1417 c------------------------------------------------------------------------------
1418       subroutine hpb_partition
1419       implicit real*8 (a-h,o-z)
1420       include 'DIMENSIONS'
1421 #ifdef MPI
1422       include 'mpif.h'
1423 #endif
1424       include 'COMMON.SBRIDGE'
1425       include 'COMMON.IOUNITS'
1426       include 'COMMON.SETUP'
1427 #ifdef MPI
1428       call int_bounds(nhpb,link_start,link_end)
1429       write (iout,*) 'Processor',fg_rank,' CG group',kolor,
1430      &  ' absolute rank',MyRank,
1431      &  ' nhpb',nhpb,' link_start=',link_start,
1432      &  ' link_end',link_end
1433 #else
1434       link_start=1
1435       link_end=nhpb
1436 #endif
1437       return
1438       end