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       ilip_end=ilip_end
 665       call int_bounds(nres-1,itube_start,itube_end)
 666       itube_start=itube_start
 667       itube_end=itube_end
 668       call int_bounds1(nres-1,ivec_start,ivec_end)
 669 c      print *,"Processor",myrank,fg_rank,fg_rank1,
 670 c     &  " ivec_start",ivec_start," ivec_end",ivec_end
 671       iset_start=loc_start+2
 672       iset_end=loc_end+2
 673       if (ndih_constr.eq.0) then
 674         idihconstr_start=1
 675         idihconstr_end=0
 676       else
 677         call int_bounds(ndih_constr,idihconstr_start,idihconstr_end)
 678       endif
 679       if (ntheta_constr.eq.0) then
 680         ithetaconstr_start=1
 681         ithetaconstr_end=0
 682       else
 683         call int_bounds
 684      &  (ntheta_constr,ithetaconstr_start,ithetaconstr_end)
 685       endif
 686 c      nsumgrad=(nres-nnt)*(nres-nnt+1)/2
 687 c      nlen=nres-nnt+1
 688       nsumgrad=(nres-nnt)*(nres-nnt+1)/2
 689       nlen=nres-nnt+1
 690       call int_bounds(nsumgrad,ngrad_start,ngrad_end)
 691       igrad_start=((2*nlen+1)
 692      &    -sqrt(float((2*nlen-1)**2-8*(ngrad_start-1))))/2
 693       jgrad_start(igrad_start)=
 694      &    ngrad_start-(2*nlen-igrad_start)*(igrad_start-1)/2
 695      &    +igrad_start
 696       jgrad_end(igrad_start)=nres
 697       igrad_end=((2*nlen+1)
 698      &    -sqrt(float((2*nlen-1)**2-8*(ngrad_end-1))))/2
 699       if (igrad_end.gt.igrad_start) jgrad_start(igrad_end)=igrad_end+1
 700       jgrad_end(igrad_end)=ngrad_end-(2*nlen-igrad_end)*(igrad_end-1)/2
 701      &    +igrad_end
 702       do i=igrad_start+1,igrad_end-1
 703         jgrad_start(i)=i+1
 704         jgrad_end(i)=nres
 705       enddo
 706       if (lprint) then
 707         write (*,*) 'Processor:',fg_rank,' CG group',kolor,
 708      & ' absolute rank',myrank,
 709      & ' loc_start',loc_start,' loc_end',loc_end,
 710      & ' ithet_start',ithet_start,' ithet_end',ithet_end,
 711      & ' iphi_start',iphi_start,' iphi_end',iphi_end,
 712      & ' iphid_start',iphid_start,' iphid_end',iphid_end,
 713      & ' ibond_start',ibond_start,' ibond_end',ibond_end,
 714      & ' ibondp_start',ibondp_start,' ibondp_end',ibondp_end,
 715      & ' iturn3_start',iturn3_start,' iturn3_end',iturn3_end,
 716      & ' iturn4_start',iturn4_start,' iturn4_end',iturn4_end,
 717      & ' ivec_start',ivec_start,' ivec_end',ivec_end,
 718      & ' iset_start',iset_start,' iset_end',iset_end,
 719      & ' idihconstr_start',idihconstr_start,' idihconstr_end',
 720      &   idihconstr_end,
 721      & ' ithetaconstr_start',ithetaconstr_start,' ithetaconstr_end',
 722      &   ithetaconstr_end
 723
 724        write (*,*) 'Processor:',fg_rank,myrank,' igrad_start',
 725      &   igrad_start,' igrad_end',igrad_end,' ngrad_start',ngrad_start,
 726      &   ' ngrad_end',ngrad_end
 727        do i=igrad_start,igrad_end
 728          write(*,*) 'Processor:',fg_rank,myrank,i,
 729      &    jgrad_start(i),jgrad_end(i)
 730        enddo
 731       endif
 732       if (nfgtasks.gt.1) then
 733         call MPI_Allgather(ivec_start,1,MPI_INTEGER,ivec_displ(0),1,
 734      &    MPI_INTEGER,FG_COMM1,IERROR)
 735         iaux=ivec_end-ivec_start+1
 736         call MPI_Allgather(iaux,1,MPI_INTEGER,ivec_count(0),1,
 737      &    MPI_INTEGER,FG_COMM1,IERROR)
 738         call MPI_Allgather(iset_start-2,1,MPI_INTEGER,iset_displ(0),1,
 739      &    MPI_INTEGER,FG_COMM,IERROR)
 740         iaux=iset_end-iset_start+1
 741         call MPI_Allgather(iaux,1,MPI_INTEGER,iset_count(0),1,
 742      &    MPI_INTEGER,FG_COMM,IERROR)
 743         call MPI_Allgather(ibond_start,1,MPI_INTEGER,ibond_displ(0),1,
 744      &    MPI_INTEGER,FG_COMM,IERROR)
 745         iaux=ibond_end-ibond_start+1
 746         call MPI_Allgather(iaux,1,MPI_INTEGER,ibond_count(0),1,
 747      &    MPI_INTEGER,FG_COMM,IERROR)
 748         call MPI_Allgather(ithet_start,1,MPI_INTEGER,ithet_displ(0),1,
 749      &    MPI_INTEGER,FG_COMM,IERROR)
 750         iaux=ithet_end-ithet_start+1
 751         call MPI_Allgather(iaux,1,MPI_INTEGER,ithet_count(0),1,
 752      &    MPI_INTEGER,FG_COMM,IERROR)
 753         call MPI_Allgather(iphi_start,1,MPI_INTEGER,iphi_displ(0),1,
 754      &    MPI_INTEGER,FG_COMM,IERROR)
 755         iaux=iphi_end-iphi_start+1
 756         call MPI_Allgather(iaux,1,MPI_INTEGER,iphi_count(0),1,
 757      &    MPI_INTEGER,FG_COMM,IERROR)
 758         call MPI_Allgather(iphi1_start,1,MPI_INTEGER,iphi1_displ(0),1,
 759      &    MPI_INTEGER,FG_COMM,IERROR)
 760         iaux=iphi1_end-iphi1_start+1
 761         call MPI_Allgather(iaux,1,MPI_INTEGER,iphi1_count(0),1,
 762      &    MPI_INTEGER,FG_COMM,IERROR)
 763         do i=0,max_fg_procs-1
 764           do j=1,maxres
 765             ielstart_all(j,i)=0
 766             ielend_all(j,i)=0
 767           enddo
 768         enddo
 769         call MPI_Allgather(iturn3_start,1,MPI_INTEGER,
 770      &    iturn3_start_all(0),1,MPI_INTEGER,FG_COMM,IERROR)
 771         call MPI_Allgather(iturn4_start,1,MPI_INTEGER,
 772      &    iturn4_start_all(0),1,MPI_INTEGER,FG_COMM,IERROR)
 773         call MPI_Allgather(iturn3_end,1,MPI_INTEGER,
 774      &    iturn3_end_all(0),1,MPI_INTEGER,FG_COMM,IERROR)
 775         call MPI_Allgather(iturn4_end,1,MPI_INTEGER,
 776      &    iturn4_end_all(0),1,MPI_INTEGER,FG_COMM,IERROR)
 777         call MPI_Allgather(iatel_s,1,MPI_INTEGER,
 778      &    iatel_s_all(0),1,MPI_INTEGER,FG_COMM,IERROR)
 779         call MPI_Allgather(iatel_e,1,MPI_INTEGER,
 780      &    iatel_e_all(0),1,MPI_INTEGER,FG_COMM,IERROR)
 781         call MPI_Allgather(ielstart(1),maxres,MPI_INTEGER,
 782      &    ielstart_all(1,0),maxres,MPI_INTEGER,FG_COMM,IERROR)
 783         call MPI_Allgather(ielend(1),maxres,MPI_INTEGER,
 784      &    ielend_all(1,0),maxres,MPI_INTEGER,FG_COMM,IERROR)
 785         if (lprint) then
 786         write (iout,*) "iatel_s_all",(iatel_s_all(i),i=0,nfgtasks)
 787         write (iout,*) "iatel_e_all",(iatel_e_all(i),i=0,nfgtasks)
 788         write (iout,*) "iturn3_start_all",
 789      &    (iturn3_start_all(i),i=0,nfgtasks-1)
 790         write (iout,*) "iturn3_end_all",
 791      &    (iturn3_end_all(i),i=0,nfgtasks-1)
 792         write (iout,*) "iturn4_start_all",
 793      &    (iturn4_start_all(i),i=0,nfgtasks-1)
 794         write (iout,*) "iturn4_end_all",
 795      &    (iturn4_end_all(i),i=0,nfgtasks-1)
 796         write (iout,*) "The ielstart_all array"
 797         do i=nnt,nct
 798           write (iout,'(20i4)') i,(ielstart_all(i,j),j=0,nfgtasks-1)
 799         enddo
 800         write (iout,*) "The ielend_all array"
 801         do i=nnt,nct
 802           write (iout,'(20i4)') i,(ielend_all(i,j),j=0,nfgtasks-1)
 803         enddo
 804         call flush(iout)
 805         endif
 806         ntask_cont_from=0
 807         ntask_cont_to=0
 808         itask_cont_from(0)=fg_rank
 809         itask_cont_to(0)=fg_rank
 810         flag=.false.
 811         do ii=iturn3_start,iturn3_end
 812           call add_int(ii,ii+2,iturn3_sent(1,ii),
 813      &                 ntask_cont_to,itask_cont_to,flag)
 814         enddo
 815         do ii=iturn4_start,iturn4_end
 816           call add_int(ii,ii+3,iturn4_sent(1,ii),
 817      &                 ntask_cont_to,itask_cont_to,flag)
 818         enddo
 819         do ii=iturn3_start,iturn3_end
 820           call add_int_from(ii,ii+2,ntask_cont_from,itask_cont_from)
 821         enddo
 822         do ii=iturn4_start,iturn4_end
 823           call add_int_from(ii,ii+3,ntask_cont_from,itask_cont_from)
 824         enddo
 825         if (lprint) then
 826         write (iout,*) "After turn3 ntask_cont_from",ntask_cont_from,
 827      &   " ntask_cont_to",ntask_cont_to
 828         write (iout,*) "itask_cont_from",
 829      &    (itask_cont_from(i),i=1,ntask_cont_from)
 830         write (iout,*) "itask_cont_to",
 831      &    (itask_cont_to(i),i=1,ntask_cont_to)
 832         call flush(iout)
 833         endif
 834 c        write (iout,*) "Loop forward"
 835 c        call flush(iout)
 836         do i=iatel_s,iatel_e
 837 c          write (iout,*) "from loop i=",i
 838 c          call flush(iout)
 839           do j=ielstart(i),ielend(i)
 840             call add_int_from(i,j,ntask_cont_from,itask_cont_from)
 841           enddo
 842         enddo
 843 c        write (iout,*) "Loop backward iatel_e-1",iatel_e-1,
 844 c     &     " iatel_e",iatel_e
 845 c        call flush(iout)
 846         nat_sent=0
 847         do i=iatel_s,iatel_e
 848 c          write (iout,*) "i",i," ielstart",ielstart(i),
 849 c     &      " ielend",ielend(i)
 850 c          call flush(iout)
 851           flag=.false.
 852           do j=ielstart(i),ielend(i)
 853             call add_int(i,j,iint_sent(1,j,nat_sent+1),ntask_cont_to,
 854      &                   itask_cont_to,flag)
 855           enddo
 856           if (flag) then
 857             nat_sent=nat_sent+1
 858             iat_sent(nat_sent)=i
 859           endif
 860         enddo
 861         if (lprint) then
 862         write (iout,*)"After longrange ntask_cont_from",ntask_cont_from,
 863      &   " ntask_cont_to",ntask_cont_to
 864         write (iout,*) "itask_cont_from",
 865      &    (itask_cont_from(i),i=1,ntask_cont_from)
 866         write (iout,*) "itask_cont_to",
 867      &    (itask_cont_to(i),i=1,ntask_cont_to)
 868         call flush(iout)
 869         write (iout,*) "iint_sent"
 870         do i=1,nat_sent
 871           ii=iat_sent(i)
 872           write (iout,'(20i4)') ii,(j,(iint_sent(k,j,i),k=1,4),
 873      &      j=ielstart(ii),ielend(ii))
 874         enddo
 875         write (iout,*) "iturn3_sent iturn3_start",iturn3_start,
 876      &    " iturn3_end",iturn3_end
 877         write (iout,'(20i4)') (i,(iturn3_sent(j,i),j=1,4),
 878      &      i=iturn3_start,iturn3_end)
 879         write (iout,*) "iturn4_sent iturn4_start",iturn4_start,
 880      &    " iturn4_end",iturn4_end
 881         write (iout,'(20i4)') (i,(iturn4_sent(j,i),j=1,4),
 882      &      i=iturn4_start,iturn4_end)
 883         call flush(iout)
 884         endif
 885         call MPI_Gather(ntask_cont_from,1,MPI_INTEGER,
 886      &   ntask_cont_from_all,1,MPI_INTEGER,king,FG_COMM,IERR)
 887 c        write (iout,*) "Gather ntask_cont_from ended"
 888 c        call flush(iout)
 889         call MPI_Gather(itask_cont_from(0),max_fg_procs,MPI_INTEGER,
 890      &   itask_cont_from_all(0,0),max_fg_procs,MPI_INTEGER,king,
 891      &   FG_COMM,IERR)
 892 c        write (iout,*) "Gather itask_cont_from ended"
 893 c        call flush(iout)
 894         call MPI_Gather(ntask_cont_to,1,MPI_INTEGER,ntask_cont_to_all,
 895      &   1,MPI_INTEGER,king,FG_COMM,IERR)
 896 c        write (iout,*) "Gather ntask_cont_to ended"
 897 c        call flush(iout)
 898         call MPI_Gather(itask_cont_to,max_fg_procs,MPI_INTEGER,
 899      &   itask_cont_to_all,max_fg_procs,MPI_INTEGER,king,FG_COMM,IERR)
 900 c        write (iout,*) "Gather itask_cont_to ended"
 901 c        call flush(iout)
 902         if (fg_rank.eq.king) then
 903           write (iout,*)"Contact receive task map (proc, #tasks, tasks)"
 904           do i=0,nfgtasks-1
 905             write (iout,'(20i4)') i,ntask_cont_from_all(i),
 906      &       (itask_cont_from_all(j,i),j=1,ntask_cont_from_all(i))
 907           enddo
 908           write (iout,*)
 909           call flush(iout)
 910           write (iout,*) "Contact send task map (proc, #tasks, tasks)"
 911           do i=0,nfgtasks-1
 912             write (iout,'(20i4)') i,ntask_cont_to_all(i),
 913      &       (itask_cont_to_all(j,i),j=1,ntask_cont_to_all(i))
 914           enddo
 915           write (iout,*)
 916           call flush(iout)
 917 C Check if every send will have a matching receive
 918           ncheck_to=0
 919           ncheck_from=0
 920           do i=0,nfgtasks-1
 921             ncheck_to=ncheck_to+ntask_cont_to_all(i)
 922             ncheck_from=ncheck_from+ntask_cont_from_all(i)
 923           enddo
 924           write (iout,*) "Control sums",ncheck_from,ncheck_to
 925           if (ncheck_from.ne.ncheck_to) then
 926             write (iout,*) "Error: #receive differs from #send."
 927             write (iout,*) "Terminating program...!"
 928             call flush(iout)
 929             flag=.false.
 930           else
 931             flag=.true.
 932             do i=0,nfgtasks-1
 933               do j=1,ntask_cont_to_all(i)
 934                 ii=itask_cont_to_all(j,i)
 935                 do k=1,ntask_cont_from_all(ii)
 936                   if (itask_cont_from_all(k,ii).eq.i) then
 937                     if(lprint)write(iout,*)"Matching send/receive",i,ii
 938                     exit
 939                   endif
 940                 enddo
 941                 if (k.eq.ntask_cont_from_all(ii)+1) then
 942                   flag=.false.
 943                   write (iout,*) "Error: send by",j," to",ii,
 944      &            " would have no matching receive"
 945                 endif
 946               enddo
 947             enddo
 948           endif
 949           if (.not.flag) then
 950             write (iout,*) "Unmatched sends; terminating program"
 951             call flush(iout)
 952           endif
 953         endif
 954         call MPI_Bcast(flag,1,MPI_LOGICAL,king,FG_COMM,IERROR)
 955 c        write (iout,*) "flag broadcast ended flag=",flag
 956 c        call flush(iout)
 957         if (.not.flag) then
 958           call MPI_Finalize(IERROR)
 959           stop "Error in INIT_INT_TABLE: unmatched send/receive."
 960         endif
 961         call MPI_Comm_group(FG_COMM,fg_group,IERR)
 962 c        write (iout,*) "MPI_Comm_group ended"
 963 c        call flush(iout)
 964         call MPI_Group_incl(fg_group,ntask_cont_from+1,
 965      &    itask_cont_from(0),CONT_FROM_GROUP,IERR)
 966         call MPI_Group_incl(fg_group,ntask_cont_to+1,itask_cont_to(0),
 967      &    CONT_TO_GROUP,IERR)
 968         do i=1,nat_sent
 969           ii=iat_sent(i)
 970           iaux=4*(ielend(ii)-ielstart(ii)+1)
 971           call MPI_Group_translate_ranks(fg_group,iaux,
 972      &      iint_sent(1,ielstart(ii),i),CONT_TO_GROUP,
 973      &      iint_sent_local(1,ielstart(ii),i),IERR )
 974 c          write (iout,*) "Ranks translated i=",i
 975 c          call flush(iout)
 976         enddo
 977         iaux=4*(iturn3_end-iturn3_start+1)
 978         call MPI_Group_translate_ranks(fg_group,iaux,
 979      &      iturn3_sent(1,iturn3_start),CONT_TO_GROUP,
 980      &      iturn3_sent_local(1,iturn3_start),IERR)
 981         iaux=4*(iturn4_end-iturn4_start+1)
 982         call MPI_Group_translate_ranks(fg_group,iaux,
 983      &      iturn4_sent(1,iturn4_start),CONT_TO_GROUP,
 984      &      iturn4_sent_local(1,iturn4_start),IERR)
 985         if (lprint) then
 986         write (iout,*) "iint_sent_local"
 987         do i=1,nat_sent
 988           ii=iat_sent(i)
 989           write (iout,'(20i4)') ii,(j,(iint_sent_local(k,j,i),k=1,4),
 990      &      j=ielstart(ii),ielend(ii))
 991           call flush(iout)
 992         enddo
 993         write (iout,*) "iturn3_sent_local iturn3_start",iturn3_start,
 994      &    " iturn3_end",iturn3_end
 995         write (iout,'(20i4)') (i,(iturn3_sent_local(j,i),j=1,4),
 996      &      i=iturn3_start,iturn3_end)
 997         write (iout,*) "iturn4_sent_local iturn4_start",iturn4_start,
 998      &    " iturn4_end",iturn4_end
 999         write (iout,'(20i4)') (i,(iturn4_sent_local(j,i),j=1,4),
1000      &      i=iturn4_start,iturn4_end)
1001         call flush(iout)
1002         endif
1003         call MPI_Group_free(fg_group,ierr)
1004         call MPI_Group_free(cont_from_group,ierr)
1005         call MPI_Group_free(cont_to_group,ierr)
1006         call MPI_Type_contiguous(3,MPI_DOUBLE_PRECISION,MPI_UYZ,IERROR)
1007         call MPI_Type_commit(MPI_UYZ,IERROR)
1008         call MPI_Type_contiguous(maxcontsshi,MPI_INTEGER,MPI_I50,IERROR)
1009         call MPI_Type_commit(MPI_I50,IERROR)
1010         call MPI_Type_contiguous(18,MPI_DOUBLE_PRECISION,MPI_UYZGRAD,
1011      &    IERROR)
1012         call MPI_Type_commit(MPI_UYZGRAD,IERROR)
1013          impishi=maxcontsshi*3
1014         call MPI_Type_contiguous(impishi,MPI_DOUBLE_PRECISION,
1015      &   MPI_SHI,IERROR)
1016         call MPI_Type_commit(MPI_SHI,IERROR)
1017         call MPI_Type_contiguous(2,MPI_DOUBLE_PRECISION,MPI_MU,IERROR)
1018         call MPI_Type_commit(MPI_MU,IERROR)
1019         call MPI_Type_contiguous(4,MPI_DOUBLE_PRECISION,MPI_MAT1,IERROR)
1020         call MPI_Type_commit(MPI_MAT1,IERROR)
1021         call MPI_Type_contiguous(8,MPI_DOUBLE_PRECISION,MPI_MAT2,IERROR)
1022         call MPI_Type_commit(MPI_MAT2,IERROR)
1023         call MPI_Type_contiguous(6,MPI_DOUBLE_PRECISION,MPI_THET,IERROR)
1024         call MPI_Type_commit(MPI_THET,IERROR)
1025         call MPI_Type_contiguous(9,MPI_DOUBLE_PRECISION,MPI_GAM,IERROR)
1026         call MPI_Type_commit(MPI_GAM,IERROR)
1027 #ifndef MATGATHER
1028 c 9/22/08 Derived types to send matrices which appear in correlation terms
1029         do i=0,nfgtasks-1
1030           if (ivec_count(i).eq.ivec_count(0)) then
1031             lentyp(i)=0
1032           else
1033             lentyp(i)=1
1034           endif
1035         enddo
1036         do ind_typ=lentyp(0),lentyp(nfgtasks-1)
1037         if (ind_typ.eq.0) then
1038           ichunk=ivec_count(0)
1039         else
1040           ichunk=ivec_count(1)
1041         endif
1042 c        do i=1,4
1043 c          blocklengths(i)=4
1044 c        enddo
1045 c        displs(1)=0
1046 c        do i=2,4
1047 c          displs(i)=displs(i-1)+blocklengths(i-1)*maxres
1048 c        enddo
1049 c        do i=1,4
1050 c          blocklengths(i)=blocklengths(i)*ichunk
1051 c        enddo
1052 c        write (iout,*) "blocklengths and displs"
1053 c        do i=1,4
1054 c          write (iout,*) i,blocklengths(i),displs(i)
1055 c        enddo
1056 c        call flush(iout)
1057 c        call MPI_Type_indexed(4,blocklengths(1),displs(1),
1058 c     &    MPI_DOUBLE_PRECISION,MPI_ROTAT1(ind_typ),IERROR)
1059 c        call MPI_Type_commit(MPI_ROTAT1(ind_typ),IERROR)
1060 c        write (iout,*) "MPI_ROTAT1",MPI_ROTAT1
1061 c        do i=1,4
1062 c          blocklengths(i)=2
1063 c        enddo
1064 c        displs(1)=0
1065 c        do i=2,4
1066 c          displs(i)=displs(i-1)+blocklengths(i-1)*maxres
1067 c        enddo
1068 c        do i=1,4
1069 c          blocklengths(i)=blocklengths(i)*ichunk
1070 c        enddo
1071 c        write (iout,*) "blocklengths and displs"
1072 c        do i=1,4
1073 c          write (iout,*) i,blocklengths(i),displs(i)
1074 c        enddo
1075 c        call flush(iout)
1076 c        call MPI_Type_indexed(4,blocklengths(1),displs(1),
1077 c     &    MPI_DOUBLE_PRECISION,MPI_ROTAT2(ind_typ),IERROR)
1078 c        call MPI_Type_commit(MPI_ROTAT2(ind_typ),IERROR)
1079 c        write (iout,*) "MPI_ROTAT2",MPI_ROTAT2
1080         do i=1,8
1081           blocklengths(i)=2
1082         enddo
1083         displs(1)=0
1084         do i=2,8
1085           displs(i)=displs(i-1)+blocklengths(i-1)*maxres
1086         enddo
1087         do i=1,15
1088           blocklengths(i)=blocklengths(i)*ichunk
1089         enddo
1090         call MPI_Type_indexed(8,blocklengths,displs,
1091      &    MPI_DOUBLE_PRECISION,MPI_PRECOMP11(ind_typ),IERROR)
1092         call MPI_Type_commit(MPI_PRECOMP11(ind_typ),IERROR)
1093         do i=1,8
1094           blocklengths(i)=4
1095         enddo
1096         displs(1)=0
1097         do i=2,8
1098           displs(i)=displs(i-1)+blocklengths(i-1)*maxres
1099         enddo
1100         do i=1,15
1101           blocklengths(i)=blocklengths(i)*ichunk
1102         enddo
1103         call MPI_Type_indexed(8,blocklengths,displs,
1104      &    MPI_DOUBLE_PRECISION,MPI_PRECOMP12(ind_typ),IERROR)
1105         call MPI_Type_commit(MPI_PRECOMP12(ind_typ),IERROR)
1106         do i=1,6
1107           blocklengths(i)=4
1108         enddo
1109         displs(1)=0
1110         do i=2,6
1111           displs(i)=displs(i-1)+blocklengths(i-1)*maxres
1112         enddo
1113         do i=1,6
1114           blocklengths(i)=blocklengths(i)*ichunk
1115         enddo
1116         call MPI_Type_indexed(6,blocklengths,displs,
1117      &    MPI_DOUBLE_PRECISION,MPI_PRECOMP22(ind_typ),IERROR)
1118         call MPI_Type_commit(MPI_PRECOMP22(ind_typ),IERROR)
1119         do i=1,2
1120           blocklengths(i)=8
1121         enddo
1122         displs(1)=0
1123         do i=2,2
1124           displs(i)=displs(i-1)+blocklengths(i-1)*maxres
1125         enddo
1126         do i=1,2
1127           blocklengths(i)=blocklengths(i)*ichunk
1128         enddo
1129         call MPI_Type_indexed(2,blocklengths,displs,
1130      &    MPI_DOUBLE_PRECISION,MPI_PRECOMP23(ind_typ),IERROR)
1131         call MPI_Type_commit(MPI_PRECOMP23(ind_typ),IERROR)
1132         do i=1,4
1133           blocklengths(i)=1
1134         enddo
1135         displs(1)=0
1136         do i=2,4
1137           displs(i)=displs(i-1)+blocklengths(i-1)*maxres
1138         enddo
1139         do i=1,4
1140           blocklengths(i)=blocklengths(i)*ichunk
1141         enddo
1142         call MPI_Type_indexed(4,blocklengths,displs,
1143      &    MPI_DOUBLE_PRECISION,MPI_ROTAT_OLD(ind_typ),IERROR)
1144         call MPI_Type_commit(MPI_ROTAT_OLD(ind_typ),IERROR)
1145         enddo
1146 #endif
1147       endif
1148       iint_start=ivec_start+1
1149       iint_end=ivec_end+1
1150       do i=0,nfgtasks-1
1151           iint_count(i)=ivec_count(i)
1152           iint_displ(i)=ivec_displ(i)
1153           ivec_displ(i)=ivec_displ(i)-1
1154           iset_displ(i)=iset_displ(i)-1
1155           ithet_displ(i)=ithet_displ(i)-1
1156           iphi_displ(i)=iphi_displ(i)-1
1157           iphi1_displ(i)=iphi1_displ(i)-1
1158           ibond_displ(i)=ibond_displ(i)-1
1159       enddo
1160       if (nfgtasks.gt.1 .and. fg_rank.eq.king
1161      &     .and. (me.eq.0 .or. .not. out1file)) then
1162         write (iout,*) "IVEC_DISPL, IVEC_COUNT, ISET_START, ISET_COUNT"
1163         do i=0,nfgtasks-1
1164           write (iout,*) i,ivec_displ(i),ivec_count(i),iset_displ(i),
1165      &      iset_count(i)
1166         enddo
1167         write (iout,*) "iphi_start",iphi_start," iphi_end",iphi_end,
1168      &    " iphi1_start",iphi1_start," iphi1_end",iphi1_end
1169         write (iout,*)"IPHI_COUNT, IPHI_DISPL, IPHI1_COUNT, IPHI1_DISPL"
1170         do i=0,nfgtasks-1
1171           write (iout,*) i,iphi_count(i),iphi_displ(i),iphi1_count(i),
1172      &      iphi1_displ(i)
1173         enddo
1174         write(iout,'(i10,a,i10,a,i10,a/a,i3,a)') n_sc_int_tot,' SC-SC ',
1175      & nele_int_tot,' electrostatic and ',nscp_int_tot,
1176      & ' SC-p interactions','were distributed among',nfgtasks,
1177      & ' fine-grain processors.'
1178       endif
1179 #else
1180       loc_start=2
1181       loc_end=nres-1
1182       ithet_start=3
1183       ithet_end=nres
1184       iturn3_start=nnt
1185       iturn3_end=nct-3
1186       iturn4_start=nnt
1187       iturn4_end=nct-4
1188       iphi_start=nnt+3
1189       iphi_end=nct
1190       iphi1_start=4
1191       iphi1_end=nres
1192       idihconstr_start=1
1193       idihconstr_end=ndih_constr
1194       ithetaconstr_start=1
1195       ithetaconstr_end=ntheta_constr
1196       iphid_start=iphi_start
1197       iphid_end=iphi_end-1
1198       itau_start=4
1199       itau_end=nres
1200       ibond_start=2
1201       ibond_end=nres-1
1202       ibondp_start=nnt
1203 C      ibondp_end=nct-1
1204       ibondp_end=nct
1205       ivec_start=1
1206       ivec_end=nres-1
1207       iset_start=3
1208       iset_end=nres+1
1209       iint_start=2
1210       iint_end=nres-1
1211       ilip_start=1
1212       ilip_end=nres
1213       itube_start=1
1214       itube_end=nres
1215
1216 #endif
1217       return
1218       end
1219 #ifdef MPI
1220 c---------------------------------------------------------------------------
1221       subroutine add_int(ii,jj,itask,ntask_cont_to,itask_cont_to,flag)
1222       implicit none
1223       include "DIMENSIONS"
1224       include "COMMON.INTERACT"
1225       include "COMMON.SETUP"
1226       include "COMMON.IOUNITS"
1227       integer ii,jj,itask(4),ntask_cont_to,
1228      &itask_cont_to(0:max_fg_procs-1)
1229       logical flag
1230       integer iturn3_start_all,iturn3_end_all,iturn4_start_all,
1231      & iturn4_end_all,iatel_s_all,iatel_e_all,ielstart_all,ielend_all
1232       common /przechowalnia/ iturn3_start_all(0:max_fg_procs),
1233      & iturn3_end_all(0:max_fg_procs),iturn4_start_all(0:max_fg_procs),
1234      & iturn4_end_all(0:max_fg_procs),iatel_s_all(0:max_fg_procs),
1235      &iatel_e_all(0:max_fg_procs),ielstart_all(maxres,0:max_fg_procs-1),
1236      & ielend_all(maxres,0:max_fg_procs-1)
1237       integer iproc,isent,k,l
1238 c Determines whether to send interaction ii,jj to other processors; a given
1239 c interaction can be sent to at most 2 processors.
1240 c Sets flag=.true. if interaction ii,jj needs to be sent to at least
1241 c one processor, otherwise flag is unchanged from the input value.
1242       isent=0
1243       itask(1)=fg_rank
1244       itask(2)=fg_rank
1245       itask(3)=fg_rank
1246       itask(4)=fg_rank
1247 c      write (iout,*) "ii",ii," jj",jj
1248 c Loop over processors to check if anybody could need interaction ii,jj
1249       do iproc=0,fg_rank-1
1250 c Check if the interaction matches any turn3 at iproc
1251         do k=iturn3_start_all(iproc),iturn3_end_all(iproc)
1252           l=k+2
1253           if (k.eq.ii-1 .and. l.eq.jj-1 .or. k.eq.ii-1 .and. l.eq.jj+1
1254      &   .or. k.eq.ii+1 .and. l.eq.jj+1 .or. k.eq.ii+1 .and. l.eq.jj-1)
1255      &    then
1256 c            write (iout,*) "turn3 to iproc",iproc," ij",ii,jj,"kl",k,l
1257 c            call flush(iout)
1258             flag=.true.
1259             if (iproc.ne.itask(1).and.iproc.ne.itask(2)
1260      &        .and.iproc.ne.itask(3).and.iproc.ne.itask(4)) then
1261               isent=isent+1
1262               itask(isent)=iproc
1263               call add_task(iproc,ntask_cont_to,itask_cont_to)
1264             endif
1265           endif
1266         enddo
1267 C Check if the interaction matches any turn4 at iproc
1268         do k=iturn4_start_all(iproc),iturn4_end_all(iproc)
1269           l=k+3
1270           if (k.eq.ii-1 .and. l.eq.jj-1 .or. k.eq.ii-1 .and. l.eq.jj+1
1271      &   .or. k.eq.ii+1 .and. l.eq.jj+1 .or. k.eq.ii+1 .and. l.eq.jj-1)
1272      &    then
1273 c            write (iout,*) "turn3 to iproc",iproc," ij",ii,jj," kl",k,l
1274 c            call flush(iout)
1275             flag=.true.
1276             if (iproc.ne.itask(1).and.iproc.ne.itask(2)
1277      &        .and.iproc.ne.itask(3).and.iproc.ne.itask(4)) then
1278               isent=isent+1
1279               itask(isent)=iproc
1280               call add_task(iproc,ntask_cont_to,itask_cont_to)
1281             endif
1282           endif
1283         enddo
1284         if (iatel_s_all(iproc).gt.0 .and. iatel_e_all(iproc).gt.0 .and.
1285      &  iatel_s_all(iproc).le.ii-1 .and. iatel_e_all(iproc).ge.ii-1)then
1286           if (ielstart_all(ii-1,iproc).le.jj-1.and.
1287      &        ielend_all(ii-1,iproc).ge.jj-1) then
1288             flag=.true.
1289             if (iproc.ne.itask(1).and.iproc.ne.itask(2)
1290      &        .and.iproc.ne.itask(3).and.iproc.ne.itask(4)) then
1291               isent=isent+1
1292               itask(isent)=iproc
1293               call add_task(iproc,ntask_cont_to,itask_cont_to)
1294             endif
1295           endif
1296           if (ielstart_all(ii-1,iproc).le.jj+1.and.
1297      &        ielend_all(ii-1,iproc).ge.jj+1) then
1298             flag=.true.
1299             if (iproc.ne.itask(1).and.iproc.ne.itask(2)
1300      &        .and.iproc.ne.itask(3).and.iproc.ne.itask(4)) then
1301               isent=isent+1
1302               itask(isent)=iproc
1303               call add_task(iproc,ntask_cont_to,itask_cont_to)
1304             endif
1305           endif
1306         endif
1307       enddo
1308       return
1309       end
1310 c---------------------------------------------------------------------------
1311       subroutine add_int_from(ii,jj,ntask_cont_from,itask_cont_from)
1312       implicit none
1313       include "DIMENSIONS"
1314       include "COMMON.INTERACT"
1315       include "COMMON.SETUP"
1316       include "COMMON.IOUNITS"
1317       integer ii,jj,itask(2),ntask_cont_from,
1318      & itask_cont_from(0:max_fg_procs-1)
1319       logical flag
1320       integer iturn3_start_all,iturn3_end_all,iturn4_start_all,
1321      & iturn4_end_all,iatel_s_all,iatel_e_all,ielstart_all,ielend_all
1322       common /przechowalnia/ iturn3_start_all(0:max_fg_procs),
1323      & iturn3_end_all(0:max_fg_procs),iturn4_start_all(0:max_fg_procs),
1324      & iturn4_end_all(0:max_fg_procs),iatel_s_all(0:max_fg_procs),
1325      &iatel_e_all(0:max_fg_procs),ielstart_all(maxres,0:max_fg_procs-1),
1326      & ielend_all(maxres,0:max_fg_procs-1)
1327       integer iproc,k,l
1328       do iproc=fg_rank+1,nfgtasks-1
1329         do k=iturn3_start_all(iproc),iturn3_end_all(iproc)
1330           l=k+2
1331           if (k.eq.ii+1 .and. l.eq.jj+1 .or. k.eq.ii+1.and.l.eq.jj-1
1332      &   .or. k.eq.ii-1 .and. l.eq.jj-1 .or. k.eq.ii-1 .and. l.eq.jj+1)
1333      &    then
1334 c            write (iout,*)"turn3 from iproc",iproc," ij",ii,jj," kl",k,l
1335             call add_task(iproc,ntask_cont_from,itask_cont_from)
1336           endif
1337         enddo
1338         do k=iturn4_start_all(iproc),iturn4_end_all(iproc)
1339           l=k+3
1340           if (k.eq.ii+1 .and. l.eq.jj+1 .or. k.eq.ii+1.and.l.eq.jj-1
1341      &   .or. k.eq.ii-1 .and. l.eq.jj-1 .or. k.eq.ii-1 .and. l.eq.jj+1)
1342      &    then
1343 c            write (iout,*)"turn4 from iproc",iproc," ij",ii,jj," kl",k,l
1344             call add_task(iproc,ntask_cont_from,itask_cont_from)
1345           endif
1346         enddo
1347         if (iatel_s_all(iproc).gt.0 .and. iatel_e_all(iproc).gt.0) then
1348           if (ii+1.ge.iatel_s_all(iproc).and.ii+1.le.iatel_e_all(iproc))
1349      &    then
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             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           endif
1359           if (ii-1.ge.iatel_s_all(iproc).and.ii-1.le.iatel_e_all(iproc))
1360      &    then
1361             if (jj-1.ge.ielstart_all(ii-1,iproc).and.
1362      &          jj-1.le.ielend_all(ii-1,iproc)) then
1363               call add_task(iproc,ntask_cont_from,itask_cont_from)
1364             endif
1365             if (jj+1.ge.ielstart_all(ii-1,iproc).and.
1366      &          jj+1.le.ielend_all(ii-1,iproc)) then
1367                call add_task(iproc,ntask_cont_from,itask_cont_from)
1368             endif
1369           endif
1370         endif
1371       enddo
1372       return
1373       end
1374 c---------------------------------------------------------------------------
1375       subroutine add_task(iproc,ntask_cont,itask_cont)
1376       implicit none
1377       include "DIMENSIONS"
1378       integer iproc,ntask_cont,itask_cont(0:max_fg_procs-1)
1379       integer ii
1380       do ii=1,ntask_cont
1381         if (itask_cont(ii).eq.iproc) return
1382       enddo
1383       ntask_cont=ntask_cont+1
1384       itask_cont(ntask_cont)=iproc
1385       return
1386       end
1387 c---------------------------------------------------------------------------
1388       subroutine int_bounds(total_ints,lower_bound,upper_bound)
1389       implicit real*8 (a-h,o-z)
1390       include 'DIMENSIONS'
1391       include 'mpif.h'
1392       include 'COMMON.SETUP'
1393       integer total_ints,lower_bound,upper_bound
1394       integer int4proc(0:max_fg_procs),sint4proc(0:max_fg_procs)
1395       nint=total_ints/nfgtasks
1396       do i=1,nfgtasks
1397         int4proc(i-1)=nint
1398       enddo
1399       nexcess=total_ints-nint*nfgtasks
1400       do i=1,nexcess
1401         int4proc(nfgtasks-i)=int4proc(nfgtasks-i)+1
1402       enddo
1403       lower_bound=0
1404       do i=0,fg_rank-1
1405         lower_bound=lower_bound+int4proc(i)
1406       enddo
1407       upper_bound=lower_bound+int4proc(fg_rank)
1408       lower_bound=lower_bound+1
1409       return
1410       end
1411 c---------------------------------------------------------------------------
1412       subroutine int_bounds1(total_ints,lower_bound,upper_bound)
1413       implicit real*8 (a-h,o-z)
1414       include 'DIMENSIONS'
1415       include 'mpif.h'
1416       include 'COMMON.SETUP'
1417       integer total_ints,lower_bound,upper_bound
1418       integer int4proc(0:max_fg_procs),sint4proc(0:max_fg_procs)
1419       nint=total_ints/nfgtasks1
1420       do i=1,nfgtasks1
1421         int4proc(i-1)=nint
1422       enddo
1423       nexcess=total_ints-nint*nfgtasks1
1424       do i=1,nexcess
1425         int4proc(nfgtasks1-i)=int4proc(nfgtasks1-i)+1
1426       enddo
1427       lower_bound=0
1428       do i=0,fg_rank1-1
1429         lower_bound=lower_bound+int4proc(i)
1430       enddo
1431       upper_bound=lower_bound+int4proc(fg_rank1)
1432       lower_bound=lower_bound+1
1433       return
1434       end
1435 c---------------------------------------------------------------------------
1436       subroutine int_partition(int_index,lower_index,upper_index,atom,
1437      & at_start,at_end,first_atom,last_atom,int_gr,jat_start,jat_end,*)
1438       implicit real*8 (a-h,o-z)
1439       include 'DIMENSIONS'
1440       include 'COMMON.IOUNITS'
1441       integer int_index,lower_index,upper_index,atom,at_start,at_end,
1442      & first_atom,last_atom,int_gr,jat_start,jat_end
1443       logical lprn
1444       lprn=.false.
1445       if (lprn) write (iout,*) 'int_index=',int_index
1446       int_index_old=int_index
1447       int_index=int_index+last_atom-first_atom+1
1448       if (lprn)
1449      &   write (iout,*) 'int_index=',int_index,
1450      &               ' int_index_old',int_index_old,
1451      &               ' lower_index=',lower_index,
1452      &               ' upper_index=',upper_index,
1453      &               ' atom=',atom,' first_atom=',first_atom,
1454      &               ' last_atom=',last_atom
1455       if (int_index.ge.lower_index) then
1456         int_gr=int_gr+1
1457         if (at_start.eq.0) then
1458           at_start=atom
1459           jat_start=first_atom-1+lower_index-int_index_old
1460         else
1461           jat_start=first_atom
1462         endif
1463         if (lprn) write (iout,*) 'jat_start',jat_start
1464         if (int_index.ge.upper_index) then
1465           at_end=atom
1466           jat_end=first_atom-1+upper_index-int_index_old
1467           return1
1468         else
1469           jat_end=last_atom
1470         endif
1471         if (lprn) write (iout,*) 'jat_end',jat_end
1472       endif
1473       return
1474       end
1475 #endif
1476 c------------------------------------------------------------------------------
1477       subroutine hpb_partition
1478       implicit real*8 (a-h,o-z)
1479       include 'DIMENSIONS'
1480 #ifdef MPI
1481       include 'mpif.h'
1482 #endif
1483       include 'COMMON.SBRIDGE'
1484       include 'COMMON.IOUNITS'
1485       include 'COMMON.SETUP'
1486 #ifdef MPI
1487       call int_bounds(nhpb,link_start,link_end)
1488       write (iout,*) 'Processor',fg_rank,' CG group',kolor,
1489      &  ' absolute rank',MyRank,
1490      &  ' nhpb',nhpb,' link_start=',link_start,
1491      &  ' link_end',link_end
1492 #else
1493       link_start=1
1494       link_end=nhpb
1495 #endif
1496       return
1497       end