source/unres/src_MD-M/initialize_p.F

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