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