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