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 Initialize correlation arrays
 249       do i=-maxtor,maxtor
 250        do k=1,2
 251         b1(k,i)=0.0
 252         b2(k,i)=0.0
 253         b1tilde(k,i)=0.0
 254 c        b2tilde(k,i)=0.0
 255         do j=1,2
 256         CC(j,k,i)=0.0
 257         Ctilde(j,k,i)=0.0
 258         DD(j,k,i)=0.0
 259         Dtilde(j,k,i)=0.0
 260         EE(j,k,i)=0.0
 261         enddo
 262        enddo
 263       enddo
 264 C
 265 C Initialize timing.
 266 C
 267       call set_timers
 268 C
 269 C Initialize variables used in minimization.
 270 C
 271 c     maxfun=5000
 272 c     maxit=2000
 273       maxfun=500
 274       maxit=200
 275       tolf=1.0D-2
 276       rtolf=5.0D-4
 277 C
 278 C Initialize the variables responsible for the mode of gradient storage.
 279 C
 280       nfl=0
 281       icg=1
 282 C
 283 C Initialize constants used to split the energy into long- and short-range
 284 C components
 285 C
 286 C      r_cut=2.0d0
 287 C      rlamb=0.3d0
 288 #ifndef SPLITELE
 289       nprint_ene=nprint_ene-1
 290 #endif
 291       return
 292       end
 293 c-------------------------------------------------------------------------
 294       block data nazwy
 295       implicit real*8 (a-h,o-z)
 296       include 'DIMENSIONS'
 297       include 'COMMON.NAMES'
 298       include 'COMMON.FFIELD'
 299       data restyp /
 300      &'DD','DAU','DAI','DDB','DSM','DPR','DLY','DAR','DHI','DAS','DGL',
 301      & 'DSG','DGN','DSN','DTH',
 302      &'DYY','DAL','DTY','DTR','DVA','DLE','DIL','DPN','MED','DCY','ZER',
 303      &'CYS','MET','PHE','ILE','LEU','VAL','TRP','TYR','ALA','GLY','THR',
 304      &'SER','GLN','ASN','GLU','ASP','HIS','ARG','LYS','PRO','SME','DBZ',
 305      &'AIB','ABU','D'/
 306       data onelet /
 307      &'z','z','z','z','z','p','k','r','h','d','e','n','q','s','t','g',
 308      &'a','y','w','v','l','i','f','m','c','x',
 309      &'C','M','F','I','L','V','W','Y','A','G','T',
 310      &'S','Q','N','E','D','H','R','K','P','z','z','z','z','X'/
 311       data potname /'LJ','LJK','BP','GB','GBV'/
 312       data ename /
 313      &   "EVDW SC-SC","EVDW2 SC-p","EES p-p","ECORR4 ","ECORR5 ",
 314      &   "ECORR6 ","EELLO ","ETURN3 ","ETURN4 ","ETURN6 ",
 315      &   "EBE bend","ESC SCloc","ETORS ","ETORSD ","EHPB ","EVDWPP ",
 316      &   "ESTR ","EVDW2_14 ","UCONST ", "      ","ESCCOR"/
 317       data wname /
 318      &   "WSC","WSCP","WELEC","WCORR","WCORR5","WCORR6","WEL_LOC",
 319      &   "WTURN3","WTURN4","WTURN6","WANG","WSCLOC","WTOR","WTORD",
 320      &   "WSTRAIN","WVDWPP","WBOND","SCAL14","     ","    ","WSCCOR"/
 321       data nprint_ene /20/
 322       data print_order/1,2,3,11,12,13,14,4,5,6,7,8,9,10,19,18,15,17,16,
 323      & 21,0/
 324       end
 325 c---------------------------------------------------------------------------
 326       subroutine init_int_table
 327       implicit real*8 (a-h,o-z)
 328       include 'DIMENSIONS'
 329 #ifdef MPI
 330       include 'mpif.h'
 331       integer blocklengths(15),displs(15)
 332 #endif
 333       include 'COMMON.CONTROL'
 334       include 'COMMON.SETUP'
 335       include 'COMMON.CHAIN'
 336       include 'COMMON.INTERACT'
 337       include 'COMMON.LOCAL'
 338       include 'COMMON.SBRIDGE'
 339       include 'COMMON.TORCNSTR'
 340       include 'COMMON.IOUNITS'
 341       include 'COMMON.DERIV'
 342       include 'COMMON.CONTACTS'
 343       common /przechowalnia/ iturn3_start_all(0:max_fg_procs),
 344      & iturn3_end_all(0:max_fg_procs),iturn4_start_all(0:max_fg_procs),
 345      & iturn4_end_all(0:max_fg_procs),iatel_s_all(0:max_fg_procs),
 346      &iatel_e_all(0:max_fg_procs),ielstart_all(maxres,0:max_fg_procs-1),
 347      & ielend_all(maxres,0:max_fg_procs-1),
 348      & ntask_cont_from_all(0:max_fg_procs-1),
 349      & itask_cont_from_all(0:max_fg_procs-1,0:max_fg_procs-1),
 350      & ntask_cont_to_all(0:max_fg_procs-1),
 351      & itask_cont_to_all(0:max_fg_procs-1,0:max_fg_procs-1)
 352       integer FG_GROUP,CONT_FROM_GROUP,CONT_TO_GROUP
 353       logical scheck,lprint,flag
 354 #ifdef MPI
 355       integer my_sc_int(0:max_fg_Procs-1),my_sc_intt(0:max_fg_Procs),
 356      & my_ele_int(0:max_fg_Procs-1),my_ele_intt(0:max_fg_Procs)
 357 C... Determine the numbers of start and end SC-SC interaction
 358 C... to deal with by current processor.
 359       do i=0,nfgtasks-1
 360         itask_cont_from(i)=fg_rank
 361         itask_cont_to(i)=fg_rank
 362       enddo
 363       lprint=.false.
 364       if (lprint)
 365      &write (iout,*) 'INIT_INT_TABLE nres=',nres,' nnt=',nnt,' nct=',nct
 366       n_sc_int_tot=(nct-nnt+1)*(nct-nnt)/2-nss
 367       call int_bounds(n_sc_int_tot,my_sc_inds,my_sc_inde)
 368       if (lprint)
 369      &  write (iout,*) 'Processor',fg_rank,' CG group',kolor,
 370      &  ' absolute rank',MyRank,
 371      &  ' n_sc_int_tot',n_sc_int_tot,' my_sc_inds=',my_sc_inds,
 372      &  ' my_sc_inde',my_sc_inde
 373       ind_sctint=0
 374       iatsc_s=0
 375       iatsc_e=0
 376 #endif
 377 c      lprint=.false.
 378       do i=1,maxres
 379         nint_gr(i)=0
 380         nscp_gr(i)=0
 381         do j=1,maxint_gr
 382           istart(i,1)=0
 383           iend(i,1)=0
 384           ielstart(i)=0
 385           ielend(i)=0
 386           iscpstart(i,1)=0
 387           iscpend(i,1)=0
 388         enddo
 389       enddo
 390       ind_scint=0
 391       ind_scint_old=0
 392 cd    write (iout,*) 'ns=',ns,' nss=',nss,' ihpb,jhpb',
 393 cd   &   (ihpb(i),jhpb(i),i=1,nss)
 394       do i=nnt,nct-1
 395         scheck=.false.
 396         do ii=1,nss
 397           if (ihpb(ii).eq.i+nres) then
 398             scheck=.true.
 399             jj=jhpb(ii)-nres
 400             goto 10
 401           endif
 402         enddo
 403    10   continue
 404 cd      write (iout,*) 'i=',i,' scheck=',scheck,' jj=',jj
 405         if (scheck) then
 406           if (jj.eq.i+1) then
 407 #ifdef MPI
 408 c            write (iout,*) 'jj=i+1'
 409             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
 410      & iatsc_s,iatsc_e,i+2,nct,nint_gr(i),istart(i,1),iend(i,1),*12)
 411 #else
 412             nint_gr(i)=1
 413             istart(i,1)=i+2
 414             iend(i,1)=nct
 415 #endif
 416           else if (jj.eq.nct) then
 417 #ifdef MPI
 418 c            write (iout,*) 'jj=nct'
 419             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
 420      &  iatsc_s,iatsc_e,i+1,nct-1,nint_gr(i),istart(i,1),iend(i,1),*12)
 421 #else
 422             nint_gr(i)=1
 423             istart(i,1)=i+1
 424             iend(i,1)=nct-1
 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,jj-1,nint_gr(i),istart(i,1),iend(i,1),*12)
 430             ii=nint_gr(i)+1
 431             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
 432      & iatsc_s,iatsc_e,jj+1,nct,nint_gr(i),istart(i,ii),iend(i,ii),*12)
 433 #else
 434             nint_gr(i)=2
 435             istart(i,1)=i+1
 436             iend(i,1)=jj-1
 437             istart(i,2)=jj+1
 438             iend(i,2)=nct
 439 #endif
 440           endif
 441         else
 442 #ifdef MPI
 443           call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
 444      &    iatsc_s,iatsc_e,i+1,nct,nint_gr(i),istart(i,1),iend(i,1),*12)
 445 #else
 446           nint_gr(i)=1
 447           istart(i,1)=i+1
 448           iend(i,1)=nct
 449           ind_scint=ind_scint+nct-i
 450 #endif
 451         endif
 452 #ifdef MPI
 453         ind_scint_old=ind_scint
 454 #endif
 455       enddo
 456    12 continue
 457 #ifndef MPI
 458       iatsc_s=nnt
 459       iatsc_e=nct-1
 460 #endif
 461 #ifdef MPI
 462       if (lprint) write (*,*) 'Processor',fg_rank,' CG Group',kolor,
 463      &   ' absolute rank',myrank,' iatsc_s=',iatsc_s,' iatsc_e=',iatsc_e
 464 #endif
 465       if (lprint) then
 466       write (iout,'(a)') 'Interaction array:'
 467       do i=iatsc_s,iatsc_e
 468         write (iout,'(i3,2(2x,2i3))')
 469      & i,(istart(i,iint),iend(i,iint),iint=1,nint_gr(i))
 470       enddo
 471       endif
 472       ispp=4
 473 #ifdef MPI
 474 C Now partition the electrostatic-interaction array
 475       npept=nct-nnt
 476       nele_int_tot=(npept-ispp)*(npept-ispp+1)/2
 477       call int_bounds(nele_int_tot,my_ele_inds,my_ele_inde)
 478       if (lprint)
 479      & write (*,*) 'Processor',fg_rank,' CG group',kolor,
 480      &  ' absolute rank',MyRank,
 481      &  ' nele_int_tot',nele_int_tot,' my_ele_inds=',my_ele_inds,
 482      &               ' my_ele_inde',my_ele_inde
 483       iatel_s=0
 484       iatel_e=0
 485       ind_eleint=0
 486       ind_eleint_old=0
 487       do i=nnt,nct-3
 488         ijunk=0
 489         call int_partition(ind_eleint,my_ele_inds,my_ele_inde,i,
 490      &    iatel_s,iatel_e,i+ispp,nct-1,ijunk,ielstart(i),ielend(i),*13)
 491       enddo ! i
 492    13 continue
 493       if (iatel_s.eq.0) iatel_s=1
 494       nele_int_tot_vdw=(npept-2)*(npept-2+1)/2
 495 c      write (iout,*) "nele_int_tot_vdw",nele_int_tot_vdw
 496       call int_bounds(nele_int_tot_vdw,my_ele_inds_vdw,my_ele_inde_vdw)
 497 c      write (iout,*) "my_ele_inds_vdw",my_ele_inds_vdw,
 498 c     & " my_ele_inde_vdw",my_ele_inde_vdw
 499       ind_eleint_vdw=0
 500       ind_eleint_vdw_old=0
 501       iatel_s_vdw=0
 502       iatel_e_vdw=0
 503       do i=nnt,nct-3
 504         ijunk=0
 505         call int_partition(ind_eleint_vdw,my_ele_inds_vdw,
 506      &    my_ele_inde_vdw,i,
 507      &    iatel_s_vdw,iatel_e_vdw,i+2,nct-1,ijunk,ielstart_vdw(i),
 508      &    ielend_vdw(i),*15)
 509 c        write (iout,*) i," ielstart_vdw",ielstart_vdw(i),
 510 c     &   " ielend_vdw",ielend_vdw(i)
 511       enddo ! i
 512       if (iatel_s_vdw.eq.0) iatel_s_vdw=1
 513    15 continue
 514 #else
 515       iatel_s=nnt
 516       iatel_e=nct-5
 517       do i=iatel_s,iatel_e
 518         ielstart(i)=i+4
 519         ielend(i)=nct-1
 520       enddo
 521       iatel_s_vdw=nnt
 522       iatel_e_vdw=nct-3
 523       do i=iatel_s_vdw,iatel_e_vdw
 524         ielstart_vdw(i)=i+2
 525         ielend_vdw(i)=nct-1
 526       enddo
 527 #endif
 528       if (lprint) then
 529         write (*,'(a)') 'Processor',fg_rank,' CG group',kolor,
 530      &  ' absolute rank',MyRank
 531         write (iout,*) 'Electrostatic interaction array:'
 532         do i=iatel_s,iatel_e
 533           write (iout,'(i3,2(2x,2i3))') i,ielstart(i),ielend(i)
 534         enddo
 535       endif ! lprint
 536 c     iscp=3
 537       iscp=2
 538 C Partition the SC-p interaction array
 539 #ifdef MPI
 540       nscp_int_tot=(npept-iscp+1)*(npept-iscp+1)
 541       call int_bounds(nscp_int_tot,my_scp_inds,my_scp_inde)
 542       if (lprint) write (iout,*) 'Processor',fg_rank,' CG group',kolor,
 543      &  ' absolute rank',myrank,
 544      &  ' nscp_int_tot',nscp_int_tot,' my_scp_inds=',my_scp_inds,
 545      &               ' my_scp_inde',my_scp_inde
 546       iatscp_s=0
 547       iatscp_e=0
 548       ind_scpint=0
 549       ind_scpint_old=0
 550       do i=nnt,nct-1
 551         if (i.lt.nnt+iscp) then
 552 cd        write (iout,*) 'i.le.nnt+iscp'
 553           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
 554      &      iatscp_s,iatscp_e,i+iscp,nct,nscp_gr(i),iscpstart(i,1),
 555      &      iscpend(i,1),*14)
 556         else if (i.gt.nct-iscp) then
 557 cd        write (iout,*) 'i.gt.nct-iscp'
 558           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
 559      &      iatscp_s,iatscp_e,nnt,i-iscp,nscp_gr(i),iscpstart(i,1),
 560      &      iscpend(i,1),*14)
 561         else
 562           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
 563      &      iatscp_s,iatscp_e,nnt,i-iscp,nscp_gr(i),iscpstart(i,1),
 564      &      iscpend(i,1),*14)
 565           ii=nscp_gr(i)+1
 566           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
 567      &      iatscp_s,iatscp_e,i+iscp,nct,nscp_gr(i),iscpstart(i,ii),
 568      &      iscpend(i,ii),*14)
 569         endif
 570       enddo ! i
 571    14 continue
 572 #else
 573       iatscp_s=nnt
 574       iatscp_e=nct-1
 575       do i=nnt,nct-1
 576         if (i.lt.nnt+iscp) then
 577           nscp_gr(i)=1
 578           iscpstart(i,1)=i+iscp
 579           iscpend(i,1)=nct
 580         elseif (i.gt.nct-iscp) then
 581           nscp_gr(i)=1
 582           iscpstart(i,1)=nnt
 583           iscpend(i,1)=i-iscp
 584         else
 585           nscp_gr(i)=2
 586           iscpstart(i,1)=nnt
 587           iscpend(i,1)=i-iscp
 588           iscpstart(i,2)=i+iscp
 589           iscpend(i,2)=nct
 590         endif
 591       enddo ! i
 592 #endif
 593       if (lprint) then
 594         write (iout,'(a)') 'SC-p interaction array:'
 595         do i=iatscp_s,iatscp_e
 596           write (iout,'(i3,2(2x,2i3))')
 597      &         i,(iscpstart(i,j),iscpend(i,j),j=1,nscp_gr(i))
 598         enddo
 599       endif ! lprint
 600 C Partition local interactions
 601 #ifdef MPI
 602       call int_bounds(nres-2,loc_start,loc_end)
 603       loc_start=loc_start+1
 604       loc_end=loc_end+1
 605       call int_bounds(nres-2,ithet_start,ithet_end)
 606       ithet_start=ithet_start+2
 607       ithet_end=ithet_end+2
 608       call int_bounds(nct-nnt-2,iturn3_start,iturn3_end)
 609       iturn3_start=iturn3_start+nnt
 610       iphi_start=iturn3_start+2
 611       iturn3_end=iturn3_end+nnt
 612       iphi_end=iturn3_end+2
 613       iturn3_start=iturn3_start-1
 614       iturn3_end=iturn3_end-1
 615       call int_bounds(nres-3,itau_start,itau_end)
 616       itau_start=itau_start+3
 617       itau_end=itau_end+3
 618       call int_bounds(nres-3,iphi1_start,iphi1_end)
 619       iphi1_start=iphi1_start+3
 620       iphi1_end=iphi1_end+3
 621       call int_bounds(nct-nnt-3,iturn4_start,iturn4_end)
 622       iturn4_start=iturn4_start+nnt
 623       iphid_start=iturn4_start+2
 624       iturn4_end=iturn4_end+nnt
 625       iphid_end=iturn4_end+2
 626       iturn4_start=iturn4_start-1
 627       iturn4_end=iturn4_end-1
 628       call int_bounds(nres-2,ibond_start,ibond_end)
 629       ibond_start=ibond_start+1
 630       ibond_end=ibond_end+1
 631       call int_bounds(nct-nnt,ibondp_start,ibondp_end)
 632       ibondp_start=ibondp_start+nnt
 633       ibondp_end=ibondp_end+nnt
 634       call int_bounds1(nres-1,ivec_start,ivec_end)
 635 c      print *,"Processor",myrank,fg_rank,fg_rank1,
 636 c     &  " ivec_start",ivec_start," ivec_end",ivec_end
 637       iset_start=loc_start+2
 638       iset_end=loc_end+2
 639       if (ndih_constr.eq.0) then
 640         idihconstr_start=1
 641         idihconstr_end=0
 642       else
 643         call int_bounds(ndih_constr,idihconstr_start,idihconstr_end)
 644       endif
 645 c      nsumgrad=(nres-nnt)*(nres-nnt+1)/2
 646 c      nlen=nres-nnt+1
 647       nsumgrad=(nres-nnt)*(nres-nnt+1)/2
 648       nlen=nres-nnt+1
 649       call int_bounds(nsumgrad,ngrad_start,ngrad_end)
 650       igrad_start=((2*nlen+1)
 651      &    -sqrt(float((2*nlen-1)**2-8*(ngrad_start-1))))/2
 652       jgrad_start(igrad_start)=
 653      &    ngrad_start-(2*nlen-igrad_start)*(igrad_start-1)/2
 654      &    +igrad_start
 655       jgrad_end(igrad_start)=nres
 656       igrad_end=((2*nlen+1)
 657      &    -sqrt(float((2*nlen-1)**2-8*(ngrad_end-1))))/2
 658       if (igrad_end.gt.igrad_start) jgrad_start(igrad_end)=igrad_end+1
 659       jgrad_end(igrad_end)=ngrad_end-(2*nlen-igrad_end)*(igrad_end-1)/2
 660      &    +igrad_end
 661       do i=igrad_start+1,igrad_end-1
 662         jgrad_start(i)=i+1
 663         jgrad_end(i)=nres
 664       enddo
 665       if (lprint) then
 666         write (*,*) 'Processor:',fg_rank,' CG group',kolor,
 667      & ' absolute rank',myrank,
 668      & ' loc_start',loc_start,' loc_end',loc_end,
 669      & ' ithet_start',ithet_start,' ithet_end',ithet_end,
 670      & ' iphi_start',iphi_start,' iphi_end',iphi_end,
 671      & ' iphid_start',iphid_start,' iphid_end',iphid_end,
 672      & ' ibond_start',ibond_start,' ibond_end',ibond_end,
 673      & ' ibondp_start',ibondp_start,' ibondp_end',ibondp_end,
 674      & ' iturn3_start',iturn3_start,' iturn3_end',iturn3_end,
 675      & ' iturn4_start',iturn4_start,' iturn4_end',iturn4_end,
 676      & ' ivec_start',ivec_start,' ivec_end',ivec_end,
 677      & ' iset_start',iset_start,' iset_end',iset_end,
 678      & ' idihconstr_start',idihconstr_start,' idihconstr_end',
 679      &   idihconstr_end
 680        write (*,*) 'Processor:',fg_rank,myrank,' igrad_start',
 681      &   igrad_start,' igrad_end',igrad_end,' ngrad_start',ngrad_start,
 682      &   ' ngrad_end',ngrad_end
 683        do i=igrad_start,igrad_end
 684          write(*,*) 'Processor:',fg_rank,myrank,i,
 685      &    jgrad_start(i),jgrad_end(i)
 686        enddo
 687       endif
 688       if (nfgtasks.gt.1) then
 689         call MPI_Allgather(ivec_start,1,MPI_INTEGER,ivec_displ(0),1,
 690      &    MPI_INTEGER,FG_COMM1,IERROR)
 691         iaux=ivec_end-ivec_start+1
 692         call MPI_Allgather(iaux,1,MPI_INTEGER,ivec_count(0),1,
 693      &    MPI_INTEGER,FG_COMM1,IERROR)
 694         call MPI_Allgather(iset_start-2,1,MPI_INTEGER,iset_displ(0),1,
 695      &    MPI_INTEGER,FG_COMM,IERROR)
 696         iaux=iset_end-iset_start+1
 697         call MPI_Allgather(iaux,1,MPI_INTEGER,iset_count(0),1,
 698      &    MPI_INTEGER,FG_COMM,IERROR)
 699         call MPI_Allgather(ibond_start,1,MPI_INTEGER,ibond_displ(0),1,
 700      &    MPI_INTEGER,FG_COMM,IERROR)
 701         iaux=ibond_end-ibond_start+1
 702         call MPI_Allgather(iaux,1,MPI_INTEGER,ibond_count(0),1,
 703      &    MPI_INTEGER,FG_COMM,IERROR)
 704         call MPI_Allgather(ithet_start,1,MPI_INTEGER,ithet_displ(0),1,
 705      &    MPI_INTEGER,FG_COMM,IERROR)
 706         iaux=ithet_end-ithet_start+1
 707         call MPI_Allgather(iaux,1,MPI_INTEGER,ithet_count(0),1,
 708      &    MPI_INTEGER,FG_COMM,IERROR)
 709         call MPI_Allgather(iphi_start,1,MPI_INTEGER,iphi_displ(0),1,
 710      &    MPI_INTEGER,FG_COMM,IERROR)
 711         iaux=iphi_end-iphi_start+1
 712         call MPI_Allgather(iaux,1,MPI_INTEGER,iphi_count(0),1,
 713      &    MPI_INTEGER,FG_COMM,IERROR)
 714         call MPI_Allgather(iphi1_start,1,MPI_INTEGER,iphi1_displ(0),1,
 715      &    MPI_INTEGER,FG_COMM,IERROR)
 716         iaux=iphi1_end-iphi1_start+1
 717         call MPI_Allgather(iaux,1,MPI_INTEGER,iphi1_count(0),1,
 718      &    MPI_INTEGER,FG_COMM,IERROR)
 719         do i=0,maxprocs-1
 720           do j=1,maxres
 721             ielstart_all(j,i)=0
 722             ielend_all(j,i)=0
 723           enddo
 724         enddo
 725         call MPI_Allgather(iturn3_start,1,MPI_INTEGER,
 726      &    iturn3_start_all(0),1,MPI_INTEGER,FG_COMM,IERROR)
 727         call MPI_Allgather(iturn4_start,1,MPI_INTEGER,
 728      &    iturn4_start_all(0),1,MPI_INTEGER,FG_COMM,IERROR)
 729         call MPI_Allgather(iturn3_end,1,MPI_INTEGER,
 730      &    iturn3_end_all(0),1,MPI_INTEGER,FG_COMM,IERROR)
 731         call MPI_Allgather(iturn4_end,1,MPI_INTEGER,
 732      &    iturn4_end_all(0),1,MPI_INTEGER,FG_COMM,IERROR)
 733         call MPI_Allgather(iatel_s,1,MPI_INTEGER,
 734      &    iatel_s_all(0),1,MPI_INTEGER,FG_COMM,IERROR)
 735         call MPI_Allgather(iatel_e,1,MPI_INTEGER,
 736      &    iatel_e_all(0),1,MPI_INTEGER,FG_COMM,IERROR)
 737         call MPI_Allgather(ielstart(1),maxres,MPI_INTEGER,
 738      &    ielstart_all(1,0),maxres,MPI_INTEGER,FG_COMM,IERROR)
 739         call MPI_Allgather(ielend(1),maxres,MPI_INTEGER,
 740      &    ielend_all(1,0),maxres,MPI_INTEGER,FG_COMM,IERROR)
 741         if (lprint) then
 742         write (iout,*) "iatel_s_all",(iatel_s_all(i),i=0,nfgtasks)
 743         write (iout,*) "iatel_e_all",(iatel_e_all(i),i=0,nfgtasks)
 744         write (iout,*) "iturn3_start_all",
 745      &    (iturn3_start_all(i),i=0,nfgtasks-1)
 746         write (iout,*) "iturn3_end_all",
 747      &    (iturn3_end_all(i),i=0,nfgtasks-1)
 748         write (iout,*) "iturn4_start_all",
 749      &    (iturn4_start_all(i),i=0,nfgtasks-1)
 750         write (iout,*) "iturn4_end_all",
 751      &    (iturn4_end_all(i),i=0,nfgtasks-1)
 752         write (iout,*) "The ielstart_all array"
 753         do i=nnt,nct
 754           write (iout,'(20i4)') i,(ielstart_all(i,j),j=0,nfgtasks-1)
 755         enddo
 756         write (iout,*) "The ielend_all array"
 757         do i=nnt,nct
 758           write (iout,'(20i4)') i,(ielend_all(i,j),j=0,nfgtasks-1)
 759         enddo
 760         call flush(iout)
 761         endif
 762         ntask_cont_from=0
 763         ntask_cont_to=0
 764         itask_cont_from(0)=fg_rank
 765         itask_cont_to(0)=fg_rank
 766         flag=.false.
 767         do ii=iturn3_start,iturn3_end
 768           call add_int(ii,ii+2,iturn3_sent(1,ii),
 769      &                 ntask_cont_to,itask_cont_to,flag)
 770         enddo
 771         do ii=iturn4_start,iturn4_end
 772           call add_int(ii,ii+3,iturn4_sent(1,ii),
 773      &                 ntask_cont_to,itask_cont_to,flag)
 774         enddo
 775         do ii=iturn3_start,iturn3_end
 776           call add_int_from(ii,ii+2,ntask_cont_from,itask_cont_from)
 777         enddo
 778         do ii=iturn4_start,iturn4_end
 779           call add_int_from(ii,ii+3,ntask_cont_from,itask_cont_from)
 780         enddo
 781         if (lprint) then
 782         write (iout,*) "After turn3 ntask_cont_from",ntask_cont_from,
 783      &   " ntask_cont_to",ntask_cont_to
 784         write (iout,*) "itask_cont_from",
 785      &    (itask_cont_from(i),i=1,ntask_cont_from)
 786         write (iout,*) "itask_cont_to",
 787      &    (itask_cont_to(i),i=1,ntask_cont_to)
 788         call flush(iout)
 789         endif
 790 c        write (iout,*) "Loop forward"
 791 c        call flush(iout)
 792         do i=iatel_s,iatel_e
 793 c          write (iout,*) "from loop i=",i
 794 c          call flush(iout)
 795           do j=ielstart(i),ielend(i)
 796             call add_int_from(i,j,ntask_cont_from,itask_cont_from)
 797           enddo
 798         enddo
 799 c        write (iout,*) "Loop backward iatel_e-1",iatel_e-1,
 800 c     &     " iatel_e",iatel_e
 801 c        call flush(iout)
 802         nat_sent=0
 803         do i=iatel_s,iatel_e
 804 c          write (iout,*) "i",i," ielstart",ielstart(i),
 805 c     &      " ielend",ielend(i)
 806 c          call flush(iout)
 807           flag=.false.
 808           do j=ielstart(i),ielend(i)
 809             call add_int(i,j,iint_sent(1,j,nat_sent+1),ntask_cont_to,
 810      &                   itask_cont_to,flag)
 811           enddo
 812           if (flag) then
 813             nat_sent=nat_sent+1
 814             iat_sent(nat_sent)=i
 815           endif
 816         enddo
 817         if (lprint) then
 818         write (iout,*)"After longrange ntask_cont_from",ntask_cont_from,
 819      &   " ntask_cont_to",ntask_cont_to
 820         write (iout,*) "itask_cont_from",
 821      &    (itask_cont_from(i),i=1,ntask_cont_from)
 822         write (iout,*) "itask_cont_to",
 823      &    (itask_cont_to(i),i=1,ntask_cont_to)
 824         call flush(iout)
 825         write (iout,*) "iint_sent"
 826         do i=1,nat_sent
 827           ii=iat_sent(i)
 828           write (iout,'(20i4)') ii,(j,(iint_sent(k,j,i),k=1,4),
 829      &      j=ielstart(ii),ielend(ii))
 830         enddo
 831         write (iout,*) "iturn3_sent iturn3_start",iturn3_start,
 832      &    " iturn3_end",iturn3_end
 833         write (iout,'(20i4)') (i,(iturn3_sent(j,i),j=1,4),
 834      &      i=iturn3_start,iturn3_end)
 835         write (iout,*) "iturn4_sent iturn4_start",iturn4_start,
 836      &    " iturn4_end",iturn4_end
 837         write (iout,'(20i4)') (i,(iturn4_sent(j,i),j=1,4),
 838      &      i=iturn4_start,iturn4_end)
 839         call flush(iout)
 840         endif
 841         call MPI_Gather(ntask_cont_from,1,MPI_INTEGER,
 842      &   ntask_cont_from_all,1,MPI_INTEGER,king,FG_COMM,IERR)
 843 c        write (iout,*) "Gather ntask_cont_from ended"
 844 c        call flush(iout)
 845         call MPI_Gather(itask_cont_from(0),max_fg_procs,MPI_INTEGER,
 846      &   itask_cont_from_all(0,0),max_fg_procs,MPI_INTEGER,king,
 847      &   FG_COMM,IERR)
 848 c        write (iout,*) "Gather itask_cont_from ended"
 849 c        call flush(iout)
 850         call MPI_Gather(ntask_cont_to,1,MPI_INTEGER,ntask_cont_to_all,
 851      &   1,MPI_INTEGER,king,FG_COMM,IERR)
 852 c        write (iout,*) "Gather ntask_cont_to ended"
 853 c        call flush(iout)
 854         call MPI_Gather(itask_cont_to,max_fg_procs,MPI_INTEGER,
 855      &   itask_cont_to_all,max_fg_procs,MPI_INTEGER,king,FG_COMM,IERR)
 856 c        write (iout,*) "Gather itask_cont_to ended"
 857 c        call flush(iout)
 858         if (fg_rank.eq.king) then
 859           write (iout,*)"Contact receive task map (proc, #tasks, tasks)"
 860           do i=0,nfgtasks-1
 861             write (iout,'(20i4)') i,ntask_cont_from_all(i),
 862      &       (itask_cont_from_all(j,i),j=1,ntask_cont_from_all(i))
 863           enddo
 864           write (iout,*)
 865           call flush(iout)
 866           write (iout,*) "Contact send task map (proc, #tasks, tasks)"
 867           do i=0,nfgtasks-1
 868             write (iout,'(20i4)') i,ntask_cont_to_all(i),
 869      &       (itask_cont_to_all(j,i),j=1,ntask_cont_to_all(i))
 870           enddo
 871           write (iout,*)
 872           call flush(iout)
 873 C Check if every send will have a matching receive
 874           ncheck_to=0
 875           ncheck_from=0
 876           do i=0,nfgtasks-1
 877             ncheck_to=ncheck_to+ntask_cont_to_all(i)
 878             ncheck_from=ncheck_from+ntask_cont_from_all(i)
 879           enddo
 880           write (iout,*) "Control sums",ncheck_from,ncheck_to
 881           if (ncheck_from.ne.ncheck_to) then
 882             write (iout,*) "Error: #receive differs from #send."
 883             write (iout,*) "Terminating program...!"
 884             call flush(iout)
 885             flag=.false.
 886           else
 887             flag=.true.
 888             do i=0,nfgtasks-1
 889               do j=1,ntask_cont_to_all(i)
 890                 ii=itask_cont_to_all(j,i)
 891                 do k=1,ntask_cont_from_all(ii)
 892                   if (itask_cont_from_all(k,ii).eq.i) then
 893                     if(lprint)write(iout,*)"Matching send/receive",i,ii
 894                     exit
 895                   endif
 896                 enddo
 897                 if (k.eq.ntask_cont_from_all(ii)+1) then
 898                   flag=.false.
 899                   write (iout,*) "Error: send by",j," to",ii,
 900      &            " would have no matching receive"
 901                 endif
 902               enddo
 903             enddo
 904           endif
 905           if (.not.flag) then
 906             write (iout,*) "Unmatched sends; terminating program"
 907             call flush(iout)
 908           endif
 909         endif
 910         call MPI_Bcast(flag,1,MPI_LOGICAL,king,FG_COMM,IERROR)
 911 c        write (iout,*) "flag broadcast ended flag=",flag
 912 c        call flush(iout)
 913         if (.not.flag) then
 914           call MPI_Finalize(IERROR)
 915           stop "Error in INIT_INT_TABLE: unmatched send/receive."
 916         endif
 917         call MPI_Comm_group(FG_COMM,fg_group,IERR)
 918 c        write (iout,*) "MPI_Comm_group ended"
 919 c        call flush(iout)
 920         call MPI_Group_incl(fg_group,ntask_cont_from+1,
 921      &    itask_cont_from(0),CONT_FROM_GROUP,IERR)
 922         call MPI_Group_incl(fg_group,ntask_cont_to+1,itask_cont_to(0),
 923      &    CONT_TO_GROUP,IERR)
 924         do i=1,nat_sent
 925           ii=iat_sent(i)
 926           iaux=4*(ielend(ii)-ielstart(ii)+1)
 927           call MPI_Group_translate_ranks(fg_group,iaux,
 928      &      iint_sent(1,ielstart(ii),i),CONT_TO_GROUP,
 929      &      iint_sent_local(1,ielstart(ii),i),IERR )
 930 c          write (iout,*) "Ranks translated i=",i
 931 c          call flush(iout)
 932         enddo
 933         iaux=4*(iturn3_end-iturn3_start+1)
 934         call MPI_Group_translate_ranks(fg_group,iaux,
 935      &      iturn3_sent(1,iturn3_start),CONT_TO_GROUP,
 936      &      iturn3_sent_local(1,iturn3_start),IERR)
 937         iaux=4*(iturn4_end-iturn4_start+1)
 938         call MPI_Group_translate_ranks(fg_group,iaux,
 939      &      iturn4_sent(1,iturn4_start),CONT_TO_GROUP,
 940      &      iturn4_sent_local(1,iturn4_start),IERR)
 941         if (lprint) then
 942         write (iout,*) "iint_sent_local"
 943         do i=1,nat_sent
 944           ii=iat_sent(i)
 945           write (iout,'(20i4)') ii,(j,(iint_sent_local(k,j,i),k=1,4),
 946      &      j=ielstart(ii),ielend(ii))
 947           call flush(iout)
 948         enddo
 949         write (iout,*) "iturn3_sent_local iturn3_start",iturn3_start,
 950      &    " iturn3_end",iturn3_end
 951         write (iout,'(20i4)') (i,(iturn3_sent_local(j,i),j=1,4),
 952      &      i=iturn3_start,iturn3_end)
 953         write (iout,*) "iturn4_sent_local iturn4_start",iturn4_start,
 954      &    " iturn4_end",iturn4_end
 955         write (iout,'(20i4)') (i,(iturn4_sent_local(j,i),j=1,4),
 956      &      i=iturn4_start,iturn4_end)
 957         call flush(iout)
 958         endif
 959         call MPI_Group_free(fg_group,ierr)
 960         call MPI_Group_free(cont_from_group,ierr)
 961         call MPI_Group_free(cont_to_group,ierr)
 962         call MPI_Type_contiguous(3,MPI_DOUBLE_PRECISION,MPI_UYZ,IERROR)
 963         call MPI_Type_commit(MPI_UYZ,IERROR)
 964         call MPI_Type_contiguous(18,MPI_DOUBLE_PRECISION,MPI_UYZGRAD,
 965      &    IERROR)
 966         call MPI_Type_commit(MPI_UYZGRAD,IERROR)
 967         call MPI_Type_contiguous(2,MPI_DOUBLE_PRECISION,MPI_MU,IERROR)
 968         call MPI_Type_commit(MPI_MU,IERROR)
 969         call MPI_Type_contiguous(4,MPI_DOUBLE_PRECISION,MPI_MAT1,IERROR)
 970         call MPI_Type_commit(MPI_MAT1,IERROR)
 971         call MPI_Type_contiguous(8,MPI_DOUBLE_PRECISION,MPI_MAT2,IERROR)
 972         call MPI_Type_commit(MPI_MAT2,IERROR)
 973         call MPI_Type_contiguous(6,MPI_DOUBLE_PRECISION,MPI_THET,IERROR)
 974         call MPI_Type_commit(MPI_THET,IERROR)
 975         call MPI_Type_contiguous(9,MPI_DOUBLE_PRECISION,MPI_GAM,IERROR)
 976         call MPI_Type_commit(MPI_GAM,IERROR)
 977 #ifndef MATGATHER
 978 c 9/22/08 Derived types to send matrices which appear in correlation terms
 979         do i=0,nfgtasks-1
 980           if (ivec_count(i).eq.ivec_count(0)) then
 981             lentyp(i)=0
 982           else
 983             lentyp(i)=1
 984           endif
 985         enddo
 986         do ind_typ=lentyp(0),lentyp(nfgtasks-1)
 987         if (ind_typ.eq.0) then
 988           ichunk=ivec_count(0)
 989         else
 990           ichunk=ivec_count(1)
 991         endif
 992 c        do i=1,4
 993 c          blocklengths(i)=4
 994 c        enddo
 995 c        displs(1)=0
 996 c        do i=2,4
 997 c          displs(i)=displs(i-1)+blocklengths(i-1)*maxres
 998 c        enddo
 999 c        do i=1,4
1000 c          blocklengths(i)=blocklengths(i)*ichunk
1001 c        enddo
1002 c        write (iout,*) "blocklengths and displs"
1003 c        do i=1,4
1004 c          write (iout,*) i,blocklengths(i),displs(i)
1005 c        enddo
1006 c        call flush(iout)
1007 c        call MPI_Type_indexed(4,blocklengths(1),displs(1),
1008 c     &    MPI_DOUBLE_PRECISION,MPI_ROTAT1(ind_typ),IERROR)
1009 c        call MPI_Type_commit(MPI_ROTAT1(ind_typ),IERROR)
1010 c        write (iout,*) "MPI_ROTAT1",MPI_ROTAT1
1011 c        do i=1,4
1012 c          blocklengths(i)=2
1013 c        enddo
1014 c        displs(1)=0
1015 c        do i=2,4
1016 c          displs(i)=displs(i-1)+blocklengths(i-1)*maxres
1017 c        enddo
1018 c        do i=1,4
1019 c          blocklengths(i)=blocklengths(i)*ichunk
1020 c        enddo
1021 c        write (iout,*) "blocklengths and displs"
1022 c        do i=1,4
1023 c          write (iout,*) i,blocklengths(i),displs(i)
1024 c        enddo
1025 c        call flush(iout)
1026 c        call MPI_Type_indexed(4,blocklengths(1),displs(1),
1027 c     &    MPI_DOUBLE_PRECISION,MPI_ROTAT2(ind_typ),IERROR)
1028 c        call MPI_Type_commit(MPI_ROTAT2(ind_typ),IERROR)
1029 c        write (iout,*) "MPI_ROTAT2",MPI_ROTAT2
1030         do i=1,8
1031           blocklengths(i)=2
1032         enddo
1033         displs(1)=0
1034         do i=2,8
1035           displs(i)=displs(i-1)+blocklengths(i-1)*maxres
1036         enddo
1037         do i=1,15
1038           blocklengths(i)=blocklengths(i)*ichunk
1039         enddo
1040         call MPI_Type_indexed(8,blocklengths,displs,
1041      &    MPI_DOUBLE_PRECISION,MPI_PRECOMP11(ind_typ),IERROR)
1042         call MPI_Type_commit(MPI_PRECOMP11(ind_typ),IERROR)
1043         do i=1,8
1044           blocklengths(i)=4
1045         enddo
1046         displs(1)=0
1047         do i=2,8
1048           displs(i)=displs(i-1)+blocklengths(i-1)*maxres
1049         enddo
1050         do i=1,15
1051           blocklengths(i)=blocklengths(i)*ichunk
1052         enddo
1053         call MPI_Type_indexed(8,blocklengths,displs,
1054      &    MPI_DOUBLE_PRECISION,MPI_PRECOMP12(ind_typ),IERROR)
1055         call MPI_Type_commit(MPI_PRECOMP12(ind_typ),IERROR)
1056         do i=1,6
1057           blocklengths(i)=4
1058         enddo
1059         displs(1)=0
1060         do i=2,6
1061           displs(i)=displs(i-1)+blocklengths(i-1)*maxres
1062         enddo
1063         do i=1,6
1064           blocklengths(i)=blocklengths(i)*ichunk
1065         enddo
1066         call MPI_Type_indexed(6,blocklengths,displs,
1067      &    MPI_DOUBLE_PRECISION,MPI_PRECOMP22(ind_typ),IERROR)
1068         call MPI_Type_commit(MPI_PRECOMP22(ind_typ),IERROR)
1069         do i=1,2
1070           blocklengths(i)=8
1071         enddo
1072         displs(1)=0
1073         do i=2,2
1074           displs(i)=displs(i-1)+blocklengths(i-1)*maxres
1075         enddo
1076         do i=1,2
1077           blocklengths(i)=blocklengths(i)*ichunk
1078         enddo
1079         call MPI_Type_indexed(2,blocklengths,displs,
1080      &    MPI_DOUBLE_PRECISION,MPI_PRECOMP23(ind_typ),IERROR)
1081         call MPI_Type_commit(MPI_PRECOMP23(ind_typ),IERROR)
1082         do i=1,4
1083           blocklengths(i)=1
1084         enddo
1085         displs(1)=0
1086         do i=2,4
1087           displs(i)=displs(i-1)+blocklengths(i-1)*maxres
1088         enddo
1089         do i=1,4
1090           blocklengths(i)=blocklengths(i)*ichunk
1091         enddo
1092         call MPI_Type_indexed(4,blocklengths,displs,
1093      &    MPI_DOUBLE_PRECISION,MPI_ROTAT_OLD(ind_typ),IERROR)
1094         call MPI_Type_commit(MPI_ROTAT_OLD(ind_typ),IERROR)
1095         enddo
1096 #endif
1097       endif
1098       iint_start=ivec_start+1
1099       iint_end=ivec_end+1
1100       do i=0,nfgtasks-1
1101           iint_count(i)=ivec_count(i)
1102           iint_displ(i)=ivec_displ(i)
1103           ivec_displ(i)=ivec_displ(i)-1
1104           iset_displ(i)=iset_displ(i)-1
1105           ithet_displ(i)=ithet_displ(i)-1
1106           iphi_displ(i)=iphi_displ(i)-1
1107           iphi1_displ(i)=iphi1_displ(i)-1
1108           ibond_displ(i)=ibond_displ(i)-1
1109       enddo
1110       if (nfgtasks.gt.1 .and. fg_rank.eq.king
1111      &     .and. (me.eq.0 .or. .not. out1file)) then
1112         write (iout,*) "IVEC_DISPL, IVEC_COUNT, ISET_START, ISET_COUNT"
1113         do i=0,nfgtasks-1
1114           write (iout,*) i,ivec_displ(i),ivec_count(i),iset_displ(i),
1115      &      iset_count(i)
1116         enddo
1117         write (iout,*) "iphi_start",iphi_start," iphi_end",iphi_end,
1118      &    " iphi1_start",iphi1_start," iphi1_end",iphi1_end
1119         write (iout,*)"IPHI_COUNT, IPHI_DISPL, IPHI1_COUNT, IPHI1_DISPL"
1120         do i=0,nfgtasks-1
1121           write (iout,*) i,iphi_count(i),iphi_displ(i),iphi1_count(i),
1122      &      iphi1_displ(i)
1123         enddo
1124         write(iout,'(i10,a,i10,a,i10,a/a,i3,a)') n_sc_int_tot,' SC-SC ',
1125      & nele_int_tot,' electrostatic and ',nscp_int_tot,
1126      & ' SC-p interactions','were distributed among',nfgtasks,
1127      & ' fine-grain processors.'
1128       endif
1129 #else
1130       loc_start=2
1131       loc_end=nres-1
1132       ithet_start=3
1133       ithet_end=nres
1134       iturn3_start=nnt
1135       iturn3_end=nct-3
1136       iturn4_start=nnt
1137       iturn4_end=nct-4
1138       iphi_start=nnt+3
1139       iphi_end=nct
1140       iphi1_start=4
1141       iphi1_end=nres
1142       idihconstr_start=1
1143       idihconstr_end=ndih_constr
1144       iphid_start=iphi_start
1145       iphid_end=iphi_end-1
1146       itau_start=4
1147       itau_end=nres
1148       ibond_start=2
1149       ibond_end=nres-1
1150       ibondp_start=nnt
1151       ibondp_end=nct-1
1152       ivec_start=1
1153       ivec_end=nres-1
1154       iset_start=3
1155       iset_end=nres+1
1156       iint_start=2
1157       iint_end=nres-1
1158 #endif
1159       return
1160       end
1161 #ifdef MPI
1162 c---------------------------------------------------------------------------
1163       subroutine add_int(ii,jj,itask,ntask_cont_to,itask_cont_to,flag)
1164       implicit none
1165       include "DIMENSIONS"
1166       include "COMMON.INTERACT"
1167       include "COMMON.SETUP"
1168       include "COMMON.IOUNITS"
1169       integer ii,jj,itask(4),ntask_cont_to,
1170      &itask_cont_to(0:max_fg_procs-1)
1171       logical flag
1172       integer iturn3_start_all,iturn3_end_all,iturn4_start_all,
1173      & iturn4_end_all,iatel_s_all,iatel_e_all,ielstart_all,ielend_all
1174       common /przechowalnia/ iturn3_start_all(0:max_fg_procs),
1175      & iturn3_end_all(0:max_fg_procs),iturn4_start_all(0:max_fg_procs),
1176      & iturn4_end_all(0:max_fg_procs),iatel_s_all(0:max_fg_procs),
1177      &iatel_e_all(0:max_fg_procs),ielstart_all(maxres,0:max_fg_procs-1),
1178      & ielend_all(maxres,0:max_fg_procs-1)
1179       integer iproc,isent,k,l
1180 c Determines whether to send interaction ii,jj to other processors; a given
1181 c interaction can be sent to at most 2 processors.
1182 c Sets flag=.true. if interaction ii,jj needs to be sent to at least
1183 c one processor, otherwise flag is unchanged from the input value.
1184       isent=0
1185       itask(1)=fg_rank
1186       itask(2)=fg_rank
1187       itask(3)=fg_rank
1188       itask(4)=fg_rank
1189 c      write (iout,*) "ii",ii," jj",jj
1190 c Loop over processors to check if anybody could need interaction ii,jj
1191       do iproc=0,fg_rank-1
1192 c Check if the interaction matches any turn3 at iproc
1193         do k=iturn3_start_all(iproc),iturn3_end_all(iproc)
1194           l=k+2
1195           if (k.eq.ii-1 .and. l.eq.jj-1 .or. k.eq.ii-1 .and. l.eq.jj+1
1196      &   .or. k.eq.ii+1 .and. l.eq.jj+1 .or. k.eq.ii+1 .and. l.eq.jj-1)
1197      &    then
1198 c            write (iout,*) "turn3 to iproc",iproc," ij",ii,jj,"kl",k,l
1199 c            call flush(iout)
1200             flag=.true.
1201             if (iproc.ne.itask(1).and.iproc.ne.itask(2)
1202      &        .and.iproc.ne.itask(3).and.iproc.ne.itask(4)) then
1203               isent=isent+1
1204               itask(isent)=iproc
1205               call add_task(iproc,ntask_cont_to,itask_cont_to)
1206             endif
1207           endif
1208         enddo
1209 C Check if the interaction matches any turn4 at iproc
1210         do k=iturn4_start_all(iproc),iturn4_end_all(iproc)
1211           l=k+3
1212           if (k.eq.ii-1 .and. l.eq.jj-1 .or. k.eq.ii-1 .and. l.eq.jj+1
1213      &   .or. k.eq.ii+1 .and. l.eq.jj+1 .or. k.eq.ii+1 .and. l.eq.jj-1)
1214      &    then
1215 c            write (iout,*) "turn3 to iproc",iproc," ij",ii,jj," kl",k,l
1216 c            call flush(iout)
1217             flag=.true.
1218             if (iproc.ne.itask(1).and.iproc.ne.itask(2)
1219      &        .and.iproc.ne.itask(3).and.iproc.ne.itask(4)) then
1220               isent=isent+1
1221               itask(isent)=iproc
1222               call add_task(iproc,ntask_cont_to,itask_cont_to)
1223             endif
1224           endif
1225         enddo
1226         if (iatel_s_all(iproc).gt.0 .and. iatel_e_all(iproc).gt.0 .and.
1227      &  iatel_s_all(iproc).le.ii-1 .and. iatel_e_all(iproc).ge.ii-1)then
1228           if (ielstart_all(ii-1,iproc).le.jj-1.and.
1229      &        ielend_all(ii-1,iproc).ge.jj-1) then
1230             flag=.true.
1231             if (iproc.ne.itask(1).and.iproc.ne.itask(2)
1232      &        .and.iproc.ne.itask(3).and.iproc.ne.itask(4)) then
1233               isent=isent+1
1234               itask(isent)=iproc
1235               call add_task(iproc,ntask_cont_to,itask_cont_to)
1236             endif
1237           endif
1238           if (ielstart_all(ii-1,iproc).le.jj+1.and.
1239      &        ielend_all(ii-1,iproc).ge.jj+1) then
1240             flag=.true.
1241             if (iproc.ne.itask(1).and.iproc.ne.itask(2)
1242      &        .and.iproc.ne.itask(3).and.iproc.ne.itask(4)) then
1243               isent=isent+1
1244               itask(isent)=iproc
1245               call add_task(iproc,ntask_cont_to,itask_cont_to)
1246             endif
1247           endif
1248         endif
1249       enddo
1250       return
1251       end
1252 c---------------------------------------------------------------------------
1253       subroutine add_int_from(ii,jj,ntask_cont_from,itask_cont_from)
1254       implicit none
1255       include "DIMENSIONS"
1256       include "COMMON.INTERACT"
1257       include "COMMON.SETUP"
1258       include "COMMON.IOUNITS"
1259       integer ii,jj,itask(2),ntask_cont_from,
1260      & itask_cont_from(0:max_fg_procs-1)
1261       logical flag
1262       integer iturn3_start_all,iturn3_end_all,iturn4_start_all,
1263      & iturn4_end_all,iatel_s_all,iatel_e_all,ielstart_all,ielend_all
1264       common /przechowalnia/ iturn3_start_all(0:max_fg_procs),
1265      & iturn3_end_all(0:max_fg_procs),iturn4_start_all(0:max_fg_procs),
1266      & iturn4_end_all(0:max_fg_procs),iatel_s_all(0:max_fg_procs),
1267      &iatel_e_all(0:max_fg_procs),ielstart_all(maxres,0:max_fg_procs-1),
1268      & ielend_all(maxres,0:max_fg_procs-1)
1269       integer iproc,k,l
1270       do iproc=fg_rank+1,nfgtasks-1
1271         do k=iturn3_start_all(iproc),iturn3_end_all(iproc)
1272           l=k+2
1273           if (k.eq.ii+1 .and. l.eq.jj+1 .or. k.eq.ii+1.and.l.eq.jj-1
1274      &   .or. k.eq.ii-1 .and. l.eq.jj-1 .or. k.eq.ii-1 .and. l.eq.jj+1)
1275      &    then
1276 c            write (iout,*)"turn3 from iproc",iproc," ij",ii,jj," kl",k,l
1277             call add_task(iproc,ntask_cont_from,itask_cont_from)
1278           endif
1279         enddo
1280         do k=iturn4_start_all(iproc),iturn4_end_all(iproc)
1281           l=k+3
1282           if (k.eq.ii+1 .and. l.eq.jj+1 .or. k.eq.ii+1.and.l.eq.jj-1
1283      &   .or. k.eq.ii-1 .and. l.eq.jj-1 .or. k.eq.ii-1 .and. l.eq.jj+1)
1284      &    then
1285 c            write (iout,*)"turn4 from iproc",iproc," ij",ii,jj," kl",k,l
1286             call add_task(iproc,ntask_cont_from,itask_cont_from)
1287           endif
1288         enddo
1289         if (iatel_s_all(iproc).gt.0 .and. iatel_e_all(iproc).gt.0) then
1290           if (ii+1.ge.iatel_s_all(iproc).and.ii+1.le.iatel_e_all(iproc))
1291      &    then
1292             if (jj+1.ge.ielstart_all(ii+1,iproc).and.
1293      &          jj+1.le.ielend_all(ii+1,iproc)) then
1294               call add_task(iproc,ntask_cont_from,itask_cont_from)
1295             endif
1296             if (jj-1.ge.ielstart_all(ii+1,iproc).and.
1297      &          jj-1.le.ielend_all(ii+1,iproc)) then
1298               call add_task(iproc,ntask_cont_from,itask_cont_from)
1299             endif
1300           endif
1301           if (ii-1.ge.iatel_s_all(iproc).and.ii-1.le.iatel_e_all(iproc))
1302      &    then
1303             if (jj-1.ge.ielstart_all(ii-1,iproc).and.
1304      &          jj-1.le.ielend_all(ii-1,iproc)) then
1305               call add_task(iproc,ntask_cont_from,itask_cont_from)
1306             endif
1307             if (jj+1.ge.ielstart_all(ii-1,iproc).and.
1308      &          jj+1.le.ielend_all(ii-1,iproc)) then
1309                call add_task(iproc,ntask_cont_from,itask_cont_from)
1310             endif
1311           endif
1312         endif
1313       enddo
1314       return
1315       end
1316 c---------------------------------------------------------------------------
1317       subroutine add_task(iproc,ntask_cont,itask_cont)
1318       implicit none
1319       include "DIMENSIONS"
1320       integer iproc,ntask_cont,itask_cont(0:max_fg_procs-1)
1321       integer ii
1322       do ii=1,ntask_cont
1323         if (itask_cont(ii).eq.iproc) return
1324       enddo
1325       ntask_cont=ntask_cont+1
1326       itask_cont(ntask_cont)=iproc
1327       return
1328       end
1329 c---------------------------------------------------------------------------
1330       subroutine int_bounds(total_ints,lower_bound,upper_bound)
1331       implicit real*8 (a-h,o-z)
1332       include 'DIMENSIONS'
1333       include 'mpif.h'
1334       include 'COMMON.SETUP'
1335       integer total_ints,lower_bound,upper_bound
1336       integer int4proc(0:max_fg_procs),sint4proc(0:max_fg_procs)
1337       nint=total_ints/nfgtasks
1338       do i=1,nfgtasks
1339         int4proc(i-1)=nint
1340       enddo
1341       nexcess=total_ints-nint*nfgtasks
1342       do i=1,nexcess
1343         int4proc(nfgtasks-i)=int4proc(nfgtasks-i)+1
1344       enddo
1345       lower_bound=0
1346       do i=0,fg_rank-1
1347         lower_bound=lower_bound+int4proc(i)
1348       enddo
1349       upper_bound=lower_bound+int4proc(fg_rank)
1350       lower_bound=lower_bound+1
1351       return
1352       end
1353 c---------------------------------------------------------------------------
1354       subroutine int_bounds1(total_ints,lower_bound,upper_bound)
1355       implicit real*8 (a-h,o-z)
1356       include 'DIMENSIONS'
1357       include 'mpif.h'
1358       include 'COMMON.SETUP'
1359       integer total_ints,lower_bound,upper_bound
1360       integer int4proc(0:max_fg_procs),sint4proc(0:max_fg_procs)
1361       nint=total_ints/nfgtasks1
1362       do i=1,nfgtasks1
1363         int4proc(i-1)=nint
1364       enddo
1365       nexcess=total_ints-nint*nfgtasks1
1366       do i=1,nexcess
1367         int4proc(nfgtasks1-i)=int4proc(nfgtasks1-i)+1
1368       enddo
1369       lower_bound=0
1370       do i=0,fg_rank1-1
1371         lower_bound=lower_bound+int4proc(i)
1372       enddo
1373       upper_bound=lower_bound+int4proc(fg_rank1)
1374       lower_bound=lower_bound+1
1375       return
1376       end
1377 c---------------------------------------------------------------------------
1378       subroutine int_partition(int_index,lower_index,upper_index,atom,
1379      & at_start,at_end,first_atom,last_atom,int_gr,jat_start,jat_end,*)
1380       implicit real*8 (a-h,o-z)
1381       include 'DIMENSIONS'
1382       include 'COMMON.IOUNITS'
1383       integer int_index,lower_index,upper_index,atom,at_start,at_end,
1384      & first_atom,last_atom,int_gr,jat_start,jat_end
1385       logical lprn
1386       lprn=.false.
1387       if (lprn) write (iout,*) 'int_index=',int_index
1388       int_index_old=int_index
1389       int_index=int_index+last_atom-first_atom+1
1390       if (lprn)
1391      &   write (iout,*) 'int_index=',int_index,
1392      &               ' int_index_old',int_index_old,
1393      &               ' lower_index=',lower_index,
1394      &               ' upper_index=',upper_index,
1395      &               ' atom=',atom,' first_atom=',first_atom,
1396      &               ' last_atom=',last_atom
1397       if (int_index.ge.lower_index) then
1398         int_gr=int_gr+1
1399         if (at_start.eq.0) then
1400           at_start=atom
1401           jat_start=first_atom-1+lower_index-int_index_old
1402         else
1403           jat_start=first_atom
1404         endif
1405         if (lprn) write (iout,*) 'jat_start',jat_start
1406         if (int_index.ge.upper_index) then
1407           at_end=atom
1408           jat_end=first_atom-1+upper_index-int_index_old
1409           return1
1410         else
1411           jat_end=last_atom
1412         endif
1413         if (lprn) write (iout,*) 'jat_end',jat_end
1414       endif
1415       return
1416       end
1417 #endif
1418 c------------------------------------------------------------------------------
1419       subroutine hpb_partition
1420       implicit real*8 (a-h,o-z)
1421       include 'DIMENSIONS'
1422 #ifdef MPI
1423       include 'mpif.h'
1424 #endif
1425       include 'COMMON.SBRIDGE'
1426       include 'COMMON.IOUNITS'
1427       include 'COMMON.SETUP'
1428 #ifdef MPI
1429       call int_bounds(nhpb,link_start,link_end)
1430       write (iout,*) 'Processor',fg_rank,' CG group',kolor,
1431      &  ' absolute rank',MyRank,
1432      &  ' nhpb',nhpb,' link_start=',link_start,
1433      &  ' link_end',link_end
1434 #else
1435       link_start=1
1436       link_end=nhpb
1437 #endif
1438       return
1439       end