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