source/unres/src_MD-M/initialize_p.F

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