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