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