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