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          if (me.eq.0 .or. .not. out1file) then
 885           write (iout,*)"Contact receive task map (proc, #tasks, tasks)"
 886           do i=0,nfgtasks-1
 887             write (iout,'(20i4)') i,ntask_cont_from_all(i),
 888      &       (itask_cont_from_all(j,i),j=1,ntask_cont_from_all(i))
 889           enddo
 890           write (iout,*)
 891           call flush(iout)
 892           write (iout,*) "Contact send task map (proc, #tasks, tasks)"
 893           do i=0,nfgtasks-1
 894             write (iout,'(20i4)') i,ntask_cont_to_all(i),
 895      &       (itask_cont_to_all(j,i),j=1,ntask_cont_to_all(i))
 896           enddo
 897           write (iout,*)
 898           call flush(iout)
 899          endif
 900 C Check if every send will have a matching receive
 901           ncheck_to=0
 902           ncheck_from=0
 903           do i=0,nfgtasks-1
 904             ncheck_to=ncheck_to+ntask_cont_to_all(i)
 905             ncheck_from=ncheck_from+ntask_cont_from_all(i)
 906           enddo
 907           if (me.eq.0 .or. .not. out1file)
 908      &     write (iout,*) "Control sums",ncheck_from,ncheck_to
 909           if (ncheck_from.ne.ncheck_to) then
 910             write (iout,*) "Error: #receive differs from #send."
 911             write (iout,*) "Terminating program...!"
 912             call flush(iout)
 913             flag=.false.
 914           else
 915             flag=.true.
 916             do i=0,nfgtasks-1
 917               do j=1,ntask_cont_to_all(i)
 918                 ii=itask_cont_to_all(j,i)
 919                 do k=1,ntask_cont_from_all(ii)
 920                   if (itask_cont_from_all(k,ii).eq.i) then
 921                     if(lprint)write(iout,*)"Matching send/receive",i,ii
 922                     exit
 923                   endif
 924                 enddo
 925                 if (k.eq.ntask_cont_from_all(ii)+1) then
 926                   flag=.false.
 927                   write (iout,*) "Error: send by",j," to",ii,
 928      &            " would have no matching receive"
 929                 endif
 930               enddo
 931             enddo
 932           endif
 933           if (.not.flag) then
 934             write (iout,*) "Unmatched sends; terminating program"
 935             call flush(iout)
 936           endif
 937         endif
 938         call MPI_Bcast(flag,1,MPI_LOGICAL,king,FG_COMM,IERROR)
 939 c        write (iout,*) "flag broadcast ended flag=",flag
 940 c        call flush(iout)
 941         if (.not.flag) then
 942           call MPI_Finalize(IERROR)
 943           stop "Error in INIT_INT_TABLE: unmatched send/receive."
 944         endif
 945         call MPI_Comm_group(FG_COMM,fg_group,IERR)
 946 c        write (iout,*) "MPI_Comm_group ended"
 947 c        call flush(iout)
 948         call MPI_Group_incl(fg_group,ntask_cont_from+1,
 949      &    itask_cont_from(0),CONT_FROM_GROUP,IERR)
 950         call MPI_Group_incl(fg_group,ntask_cont_to+1,itask_cont_to(0),
 951      &    CONT_TO_GROUP,IERR)
 952         do i=1,nat_sent
 953           ii=iat_sent(i)
 954           iaux=4*(ielend(ii)-ielstart(ii)+1)
 955           call MPI_Group_translate_ranks(fg_group,iaux,
 956      &      iint_sent(1,ielstart(ii),i),CONT_TO_GROUP,
 957      &      iint_sent_local(1,ielstart(ii),i),IERR )
 958 c          write (iout,*) "Ranks translated i=",i
 959 c          call flush(iout)
 960         enddo
 961         iaux=4*(iturn3_end-iturn3_start+1)
 962         call MPI_Group_translate_ranks(fg_group,iaux,
 963      &      iturn3_sent(1,iturn3_start),CONT_TO_GROUP,
 964      &      iturn3_sent_local(1,iturn3_start),IERR)
 965         iaux=4*(iturn4_end-iturn4_start+1)
 966         call MPI_Group_translate_ranks(fg_group,iaux,
 967      &      iturn4_sent(1,iturn4_start),CONT_TO_GROUP,
 968      &      iturn4_sent_local(1,iturn4_start),IERR)
 969         if (lprint) then
 970         write (iout,*) "iint_sent_local"
 971         do i=1,nat_sent
 972           ii=iat_sent(i)
 973           write (iout,'(20i4)') ii,(j,(iint_sent_local(k,j,i),k=1,4),
 974      &      j=ielstart(ii),ielend(ii))
 975           call flush(iout)
 976         enddo
 977         write (iout,*) "iturn3_sent_local iturn3_start",iturn3_start,
 978      &    " iturn3_end",iturn3_end
 979         write (iout,'(20i4)') (i,(iturn3_sent_local(j,i),j=1,4),
 980      &      i=iturn3_start,iturn3_end)
 981         write (iout,*) "iturn4_sent_local iturn4_start",iturn4_start,
 982      &    " iturn4_end",iturn4_end
 983         write (iout,'(20i4)') (i,(iturn4_sent_local(j,i),j=1,4),
 984      &      i=iturn4_start,iturn4_end)
 985         call flush(iout)
 986         endif
 987         call MPI_Group_free(fg_group,ierr)
 988         call MPI_Group_free(cont_from_group,ierr)
 989         call MPI_Group_free(cont_to_group,ierr)
 990         call MPI_Type_contiguous(3,MPI_DOUBLE_PRECISION,MPI_UYZ,IERROR)
 991         call MPI_Type_commit(MPI_UYZ,IERROR)
 992         call MPI_Type_contiguous(18,MPI_DOUBLE_PRECISION,MPI_UYZGRAD,
 993      &    IERROR)
 994         call MPI_Type_commit(MPI_UYZGRAD,IERROR)
 995         call MPI_Type_contiguous(2,MPI_DOUBLE_PRECISION,MPI_MU,IERROR)
 996         call MPI_Type_commit(MPI_MU,IERROR)
 997         call MPI_Type_contiguous(4,MPI_DOUBLE_PRECISION,MPI_MAT1,IERROR)
 998         call MPI_Type_commit(MPI_MAT1,IERROR)
 999         call MPI_Type_contiguous(8,MPI_DOUBLE_PRECISION,MPI_MAT2,IERROR)
1000         call MPI_Type_commit(MPI_MAT2,IERROR)
1001         call MPI_Type_contiguous(6,MPI_DOUBLE_PRECISION,MPI_THET,IERROR)
1002         call MPI_Type_commit(MPI_THET,IERROR)
1003         call MPI_Type_contiguous(9,MPI_DOUBLE_PRECISION,MPI_GAM,IERROR)
1004         call MPI_Type_commit(MPI_GAM,IERROR)
1005 #ifndef MATGATHER
1006 c 9/22/08 Derived types to send matrices which appear in correlation terms
1007         do i=0,nfgtasks-1
1008           if (ivec_count(i).eq.ivec_count(0)) then
1009             lentyp(i)=0
1010           else
1011             lentyp(i)=1
1012           endif
1013         enddo
1014         do ind_typ=lentyp(0),lentyp(nfgtasks-1)
1015         if (ind_typ.eq.0) then
1016           ichunk=ivec_count(0)
1017         else
1018           ichunk=ivec_count(1)
1019         endif
1020 c        do i=1,4
1021 c          blocklengths(i)=4
1022 c        enddo
1023 c        displs(1)=0
1024 c        do i=2,4
1025 c          displs(i)=displs(i-1)+blocklengths(i-1)*maxres
1026 c        enddo
1027 c        do i=1,4
1028 c          blocklengths(i)=blocklengths(i)*ichunk
1029 c        enddo
1030 c        write (iout,*) "blocklengths and displs"
1031 c        do i=1,4
1032 c          write (iout,*) i,blocklengths(i),displs(i)
1033 c        enddo
1034 c        call flush(iout)
1035 c        call MPI_Type_indexed(4,blocklengths(1),displs(1),
1036 c     &    MPI_DOUBLE_PRECISION,MPI_ROTAT1(ind_typ),IERROR)
1037 c        call MPI_Type_commit(MPI_ROTAT1(ind_typ),IERROR)
1038 c        write (iout,*) "MPI_ROTAT1",MPI_ROTAT1
1039 c        do i=1,4
1040 c          blocklengths(i)=2
1041 c        enddo
1042 c        displs(1)=0
1043 c        do i=2,4
1044 c          displs(i)=displs(i-1)+blocklengths(i-1)*maxres
1045 c        enddo
1046 c        do i=1,4
1047 c          blocklengths(i)=blocklengths(i)*ichunk
1048 c        enddo
1049 c        write (iout,*) "blocklengths and displs"
1050 c        do i=1,4
1051 c          write (iout,*) i,blocklengths(i),displs(i)
1052 c        enddo
1053 c        call flush(iout)
1054 c        call MPI_Type_indexed(4,blocklengths(1),displs(1),
1055 c     &    MPI_DOUBLE_PRECISION,MPI_ROTAT2(ind_typ),IERROR)
1056 c        call MPI_Type_commit(MPI_ROTAT2(ind_typ),IERROR)
1057 c        write (iout,*) "MPI_ROTAT2",MPI_ROTAT2
1058         do i=1,8
1059           blocklengths(i)=2
1060         enddo
1061         displs(1)=0
1062         do i=2,8
1063           displs(i)=displs(i-1)+blocklengths(i-1)*maxres
1064         enddo
1065         do i=1,15
1066           blocklengths(i)=blocklengths(i)*ichunk
1067         enddo
1068         call MPI_Type_indexed(8,blocklengths,displs,
1069      &    MPI_DOUBLE_PRECISION,MPI_PRECOMP11(ind_typ),IERROR)
1070         call MPI_Type_commit(MPI_PRECOMP11(ind_typ),IERROR)
1071         do i=1,8
1072           blocklengths(i)=4
1073         enddo
1074         displs(1)=0
1075         do i=2,8
1076           displs(i)=displs(i-1)+blocklengths(i-1)*maxres
1077         enddo
1078         do i=1,15
1079           blocklengths(i)=blocklengths(i)*ichunk
1080         enddo
1081         call MPI_Type_indexed(8,blocklengths,displs,
1082      &    MPI_DOUBLE_PRECISION,MPI_PRECOMP12(ind_typ),IERROR)
1083         call MPI_Type_commit(MPI_PRECOMP12(ind_typ),IERROR)
1084         do i=1,6
1085           blocklengths(i)=4
1086         enddo
1087         displs(1)=0
1088         do i=2,6
1089           displs(i)=displs(i-1)+blocklengths(i-1)*maxres
1090         enddo
1091         do i=1,6
1092           blocklengths(i)=blocklengths(i)*ichunk
1093         enddo
1094         call MPI_Type_indexed(6,blocklengths,displs,
1095      &    MPI_DOUBLE_PRECISION,MPI_PRECOMP22(ind_typ),IERROR)
1096         call MPI_Type_commit(MPI_PRECOMP22(ind_typ),IERROR)
1097         do i=1,2
1098           blocklengths(i)=8
1099         enddo
1100         displs(1)=0
1101         do i=2,2
1102           displs(i)=displs(i-1)+blocklengths(i-1)*maxres
1103         enddo
1104         do i=1,2
1105           blocklengths(i)=blocklengths(i)*ichunk
1106         enddo
1107         call MPI_Type_indexed(2,blocklengths,displs,
1108      &    MPI_DOUBLE_PRECISION,MPI_PRECOMP23(ind_typ),IERROR)
1109         call MPI_Type_commit(MPI_PRECOMP23(ind_typ),IERROR)
1110         do i=1,4
1111           blocklengths(i)=1
1112         enddo
1113         displs(1)=0
1114         do i=2,4
1115           displs(i)=displs(i-1)+blocklengths(i-1)*maxres
1116         enddo
1117         do i=1,4
1118           blocklengths(i)=blocklengths(i)*ichunk
1119         enddo
1120         call MPI_Type_indexed(4,blocklengths,displs,
1121      &    MPI_DOUBLE_PRECISION,MPI_ROTAT_OLD(ind_typ),IERROR)
1122         call MPI_Type_commit(MPI_ROTAT_OLD(ind_typ),IERROR)
1123         enddo
1124 #endif
1125       endif
1126       iint_start=ivec_start+1
1127       iint_end=ivec_end+1
1128       do i=0,nfgtasks-1
1129           iint_count(i)=ivec_count(i)
1130           iint_displ(i)=ivec_displ(i)
1131           ivec_displ(i)=ivec_displ(i)-1
1132           iset_displ(i)=iset_displ(i)-1
1133           ithet_displ(i)=ithet_displ(i)-1
1134           iphi_displ(i)=iphi_displ(i)-1
1135           iphi1_displ(i)=iphi1_displ(i)-1
1136           ibond_displ(i)=ibond_displ(i)-1
1137       enddo
1138       if (nfgtasks.gt.1 .and. fg_rank.eq.king
1139      &     .and. (me.eq.0 .or. .not. out1file)) then
1140         write (iout,*) "IVEC_DISPL, IVEC_COUNT, ISET_START, ISET_COUNT"
1141         do i=0,nfgtasks-1
1142           write (iout,*) i,ivec_displ(i),ivec_count(i),iset_displ(i),
1143      &      iset_count(i)
1144         enddo
1145         write (iout,*) "iphi_start",iphi_start," iphi_end",iphi_end,
1146      &    " iphi1_start",iphi1_start," iphi1_end",iphi1_end
1147         write (iout,*)"IPHI_COUNT, IPHI_DISPL, IPHI1_COUNT, IPHI1_DISPL"
1148         do i=0,nfgtasks-1
1149           write (iout,*) i,iphi_count(i),iphi_displ(i),iphi1_count(i),
1150      &      iphi1_displ(i)
1151         enddo
1152         write(iout,'(i10,a,i10,a,i10,a/a,i3,a)') n_sc_int_tot,' SC-SC ',
1153      & nele_int_tot,' electrostatic and ',nscp_int_tot,
1154      & ' SC-p interactions','were distributed among',nfgtasks,
1155      & ' fine-grain processors.'
1156       endif
1157 #else
1158       loc_start=2
1159       loc_end=nres-1
1160       ithet_start=3
1161       ithet_end=nres
1162       iturn3_start=nnt
1163       iturn3_end=nct-3
1164       iturn4_start=nnt
1165       iturn4_end=nct-4
1166       iphi_start=nnt+3
1167       iphi_end=nct
1168       iphi1_start=4
1169       iphi1_end=nres
1170       idihconstr_start=1
1171       idihconstr_end=ndih_constr
1172       iphid_start=iphi_start
1173       iphid_end=iphi_end-1
1174       itau_start=4
1175       itau_end=nres
1176       ibond_start=2
1177       ibond_end=nres-1
1178       ibondp_start=nnt
1179 C      ibondp_end=nct-1
1180       ibondp_end=nct
1181       ivec_start=1
1182       ivec_end=nres-1
1183       iset_start=3
1184       iset_end=nres+1
1185       iint_start=2
1186       iint_end=nres-1
1187       ilip_start=1
1188       ilip_end=nres
1189 #endif
1190       return
1191       end
1192 #ifdef MPI
1193 c---------------------------------------------------------------------------
1194       subroutine add_int(ii,jj,itask,ntask_cont_to,itask_cont_to,flag)
1195       implicit none
1196       include "DIMENSIONS"
1197       include "COMMON.INTERACT"
1198       include "COMMON.SETUP"
1199       include "COMMON.IOUNITS"
1200       integer ii,jj,itask(4),ntask_cont_to,
1201      &itask_cont_to(0:max_fg_procs-1)
1202       logical flag
1203       integer iturn3_start_all,iturn3_end_all,iturn4_start_all,
1204      & iturn4_end_all,iatel_s_all,iatel_e_all,ielstart_all,ielend_all
1205       common /przechowalnia/ iturn3_start_all(0:max_fg_procs),
1206      & iturn3_end_all(0:max_fg_procs),iturn4_start_all(0:max_fg_procs),
1207      & iturn4_end_all(0:max_fg_procs),iatel_s_all(0:max_fg_procs),
1208      &iatel_e_all(0:max_fg_procs),ielstart_all(maxres,0:max_fg_procs-1),
1209      & ielend_all(maxres,0:max_fg_procs-1)
1210       integer iproc,isent,k,l
1211 c Determines whether to send interaction ii,jj to other processors; a given
1212 c interaction can be sent to at most 2 processors.
1213 c Sets flag=.true. if interaction ii,jj needs to be sent to at least
1214 c one processor, otherwise flag is unchanged from the input value.
1215       isent=0
1216       itask(1)=fg_rank
1217       itask(2)=fg_rank
1218       itask(3)=fg_rank
1219       itask(4)=fg_rank
1220 c      write (iout,*) "ii",ii," jj",jj
1221 c Loop over processors to check if anybody could need interaction ii,jj
1222       do iproc=0,fg_rank-1
1223 c Check if the interaction matches any turn3 at iproc
1224         do k=iturn3_start_all(iproc),iturn3_end_all(iproc)
1225           l=k+2
1226           if (k.eq.ii-1 .and. l.eq.jj-1 .or. k.eq.ii-1 .and. l.eq.jj+1
1227      &   .or. k.eq.ii+1 .and. l.eq.jj+1 .or. k.eq.ii+1 .and. l.eq.jj-1)
1228      &    then
1229 c            write (iout,*) "turn3 to iproc",iproc," ij",ii,jj,"kl",k,l
1230 c            call flush(iout)
1231             flag=.true.
1232             if (iproc.ne.itask(1).and.iproc.ne.itask(2)
1233      &        .and.iproc.ne.itask(3).and.iproc.ne.itask(4)) then
1234               isent=isent+1
1235               itask(isent)=iproc
1236               call add_task(iproc,ntask_cont_to,itask_cont_to)
1237             endif
1238           endif
1239         enddo
1240 C Check if the interaction matches any turn4 at iproc
1241         do k=iturn4_start_all(iproc),iturn4_end_all(iproc)
1242           l=k+3
1243           if (k.eq.ii-1 .and. l.eq.jj-1 .or. k.eq.ii-1 .and. l.eq.jj+1
1244      &   .or. k.eq.ii+1 .and. l.eq.jj+1 .or. k.eq.ii+1 .and. l.eq.jj-1)
1245      &    then
1246 c            write (iout,*) "turn3 to iproc",iproc," ij",ii,jj," kl",k,l
1247 c            call flush(iout)
1248             flag=.true.
1249             if (iproc.ne.itask(1).and.iproc.ne.itask(2)
1250      &        .and.iproc.ne.itask(3).and.iproc.ne.itask(4)) then
1251               isent=isent+1
1252               itask(isent)=iproc
1253               call add_task(iproc,ntask_cont_to,itask_cont_to)
1254             endif
1255           endif
1256         enddo
1257         if (iatel_s_all(iproc).gt.0 .and. iatel_e_all(iproc).gt.0 .and.
1258      &  iatel_s_all(iproc).le.ii-1 .and. iatel_e_all(iproc).ge.ii-1)then
1259           if (ielstart_all(ii-1,iproc).le.jj-1.and.
1260      &        ielend_all(ii-1,iproc).ge.jj-1) then
1261             flag=.true.
1262             if (iproc.ne.itask(1).and.iproc.ne.itask(2)
1263      &        .and.iproc.ne.itask(3).and.iproc.ne.itask(4)) then
1264               isent=isent+1
1265               itask(isent)=iproc
1266               call add_task(iproc,ntask_cont_to,itask_cont_to)
1267             endif
1268           endif
1269           if (ielstart_all(ii-1,iproc).le.jj+1.and.
1270      &        ielend_all(ii-1,iproc).ge.jj+1) then
1271             flag=.true.
1272             if (iproc.ne.itask(1).and.iproc.ne.itask(2)
1273      &        .and.iproc.ne.itask(3).and.iproc.ne.itask(4)) then
1274               isent=isent+1
1275               itask(isent)=iproc
1276               call add_task(iproc,ntask_cont_to,itask_cont_to)
1277             endif
1278           endif
1279         endif
1280       enddo
1281       return
1282       end
1283 c---------------------------------------------------------------------------
1284       subroutine add_int_from(ii,jj,ntask_cont_from,itask_cont_from)
1285       implicit none
1286       include "DIMENSIONS"
1287       include "COMMON.INTERACT"
1288       include "COMMON.SETUP"
1289       include "COMMON.IOUNITS"
1290       integer ii,jj,itask(2),ntask_cont_from,
1291      & itask_cont_from(0:max_fg_procs-1)
1292       logical flag
1293       integer iturn3_start_all,iturn3_end_all,iturn4_start_all,
1294      & iturn4_end_all,iatel_s_all,iatel_e_all,ielstart_all,ielend_all
1295       common /przechowalnia/ iturn3_start_all(0:max_fg_procs),
1296      & iturn3_end_all(0:max_fg_procs),iturn4_start_all(0:max_fg_procs),
1297      & iturn4_end_all(0:max_fg_procs),iatel_s_all(0:max_fg_procs),
1298      &iatel_e_all(0:max_fg_procs),ielstart_all(maxres,0:max_fg_procs-1),
1299      & ielend_all(maxres,0:max_fg_procs-1)
1300       integer iproc,k,l
1301       do iproc=fg_rank+1,nfgtasks-1
1302         do k=iturn3_start_all(iproc),iturn3_end_all(iproc)
1303           l=k+2
1304           if (k.eq.ii+1 .and. l.eq.jj+1 .or. k.eq.ii+1.and.l.eq.jj-1
1305      &   .or. k.eq.ii-1 .and. l.eq.jj-1 .or. k.eq.ii-1 .and. l.eq.jj+1)
1306      &    then
1307 c            write (iout,*)"turn3 from iproc",iproc," ij",ii,jj," kl",k,l
1308             call add_task(iproc,ntask_cont_from,itask_cont_from)
1309           endif
1310         enddo
1311         do k=iturn4_start_all(iproc),iturn4_end_all(iproc)
1312           l=k+3
1313           if (k.eq.ii+1 .and. l.eq.jj+1 .or. k.eq.ii+1.and.l.eq.jj-1
1314      &   .or. k.eq.ii-1 .and. l.eq.jj-1 .or. k.eq.ii-1 .and. l.eq.jj+1)
1315      &    then
1316 c            write (iout,*)"turn4 from iproc",iproc," ij",ii,jj," kl",k,l
1317             call add_task(iproc,ntask_cont_from,itask_cont_from)
1318           endif
1319         enddo
1320         if (iatel_s_all(iproc).gt.0 .and. iatel_e_all(iproc).gt.0) then
1321           if (ii+1.ge.iatel_s_all(iproc).and.ii+1.le.iatel_e_all(iproc))
1322      &    then
1323             if (jj+1.ge.ielstart_all(ii+1,iproc).and.
1324      &          jj+1.le.ielend_all(ii+1,iproc)) then
1325               call add_task(iproc,ntask_cont_from,itask_cont_from)
1326             endif
1327             if (jj-1.ge.ielstart_all(ii+1,iproc).and.
1328      &          jj-1.le.ielend_all(ii+1,iproc)) then
1329               call add_task(iproc,ntask_cont_from,itask_cont_from)
1330             endif
1331           endif
1332           if (ii-1.ge.iatel_s_all(iproc).and.ii-1.le.iatel_e_all(iproc))
1333      &    then
1334             if (jj-1.ge.ielstart_all(ii-1,iproc).and.
1335      &          jj-1.le.ielend_all(ii-1,iproc)) then
1336               call add_task(iproc,ntask_cont_from,itask_cont_from)
1337             endif
1338             if (jj+1.ge.ielstart_all(ii-1,iproc).and.
1339      &          jj+1.le.ielend_all(ii-1,iproc)) then
1340                call add_task(iproc,ntask_cont_from,itask_cont_from)
1341             endif
1342           endif
1343         endif
1344       enddo
1345       return
1346       end
1347 c---------------------------------------------------------------------------
1348       subroutine add_task(iproc,ntask_cont,itask_cont)
1349       implicit none
1350       include "DIMENSIONS"
1351       integer iproc,ntask_cont,itask_cont(0:max_fg_procs-1)
1352       integer ii
1353       do ii=1,ntask_cont
1354         if (itask_cont(ii).eq.iproc) return
1355       enddo
1356       ntask_cont=ntask_cont+1
1357       itask_cont(ntask_cont)=iproc
1358       return
1359       end
1360 c---------------------------------------------------------------------------
1361       subroutine int_bounds(total_ints,lower_bound,upper_bound)
1362       implicit real*8 (a-h,o-z)
1363       include 'DIMENSIONS'
1364       include 'mpif.h'
1365       include 'COMMON.SETUP'
1366       integer total_ints,lower_bound,upper_bound
1367       integer int4proc(0:max_fg_procs),sint4proc(0:max_fg_procs)
1368       nint=total_ints/nfgtasks
1369       do i=1,nfgtasks
1370         int4proc(i-1)=nint
1371       enddo
1372       nexcess=total_ints-nint*nfgtasks
1373       do i=1,nexcess
1374         int4proc(nfgtasks-i)=int4proc(nfgtasks-i)+1
1375       enddo
1376       lower_bound=0
1377       do i=0,fg_rank-1
1378         lower_bound=lower_bound+int4proc(i)
1379       enddo
1380       upper_bound=lower_bound+int4proc(fg_rank)
1381       lower_bound=lower_bound+1
1382       return
1383       end
1384 c---------------------------------------------------------------------------
1385       subroutine int_bounds1(total_ints,lower_bound,upper_bound)
1386       implicit real*8 (a-h,o-z)
1387       include 'DIMENSIONS'
1388       include 'mpif.h'
1389       include 'COMMON.SETUP'
1390       integer total_ints,lower_bound,upper_bound
1391       integer int4proc(0:max_fg_procs),sint4proc(0:max_fg_procs)
1392       nint=total_ints/nfgtasks1
1393       do i=1,nfgtasks1
1394         int4proc(i-1)=nint
1395       enddo
1396       nexcess=total_ints-nint*nfgtasks1
1397       do i=1,nexcess
1398         int4proc(nfgtasks1-i)=int4proc(nfgtasks1-i)+1
1399       enddo
1400       lower_bound=0
1401       do i=0,fg_rank1-1
1402         lower_bound=lower_bound+int4proc(i)
1403       enddo
1404       upper_bound=lower_bound+int4proc(fg_rank1)
1405       lower_bound=lower_bound+1
1406       return
1407       end
1408 c---------------------------------------------------------------------------
1409       subroutine int_partition(int_index,lower_index,upper_index,atom,
1410      & at_start,at_end,first_atom,last_atom,int_gr,jat_start,jat_end,*)
1411       implicit real*8 (a-h,o-z)
1412       include 'DIMENSIONS'
1413       include 'COMMON.IOUNITS'
1414       integer int_index,lower_index,upper_index,atom,at_start,at_end,
1415      & first_atom,last_atom,int_gr,jat_start,jat_end
1416       logical lprn
1417       lprn=.false.
1418       if (lprn) write (iout,*) 'int_index=',int_index
1419       int_index_old=int_index
1420       int_index=int_index+last_atom-first_atom+1
1421       if (lprn)
1422      &   write (iout,*) 'int_index=',int_index,
1423      &               ' int_index_old',int_index_old,
1424      &               ' lower_index=',lower_index,
1425      &               ' upper_index=',upper_index,
1426      &               ' atom=',atom,' first_atom=',first_atom,
1427      &               ' last_atom=',last_atom
1428       if (int_index.ge.lower_index) then
1429         int_gr=int_gr+1
1430         if (at_start.eq.0) then
1431           at_start=atom
1432           jat_start=first_atom-1+lower_index-int_index_old
1433         else
1434           jat_start=first_atom
1435         endif
1436         if (lprn) write (iout,*) 'jat_start',jat_start
1437         if (int_index.ge.upper_index) then
1438           at_end=atom
1439           jat_end=first_atom-1+upper_index-int_index_old
1440           return1
1441         else
1442           jat_end=last_atom
1443         endif
1444         if (lprn) write (iout,*) 'jat_end',jat_end
1445       endif
1446       return
1447       end
1448 #endif
1449 c------------------------------------------------------------------------------
1450       subroutine hpb_partition
1451       implicit real*8 (a-h,o-z)
1452       include 'DIMENSIONS'
1453 #ifdef MPI
1454       include 'mpif.h'
1455 #endif
1456       include 'COMMON.SBRIDGE'
1457       include 'COMMON.IOUNITS'
1458       include 'COMMON.SETUP'
1459       include 'COMMON.CONTROL'
1460 c      write(2,*)"hpb_partition: nhpb=",nhpb
1461 #ifdef MPI
1462       call int_bounds(nhpb,link_start,link_end)
1463       if (.not. out1file)
1464      &  write (iout,*) 'Processor',fg_rank,' CG group',kolor,
1465      &  ' absolute rank',MyRank,
1466      &  ' nhpb',nhpb,' link_start=',link_start,
1467      &  ' link_end',link_end
1468 #else
1469       link_start=1
1470       link_end=nhpb
1471 #endif
1472 c      write(2,*)"hpb_partition: link_start=",nhpb," link_end=",link_end
1473       return
1474       end
1475 c------------------------------------------------------------------------------
1476       subroutine homology_partition
1477       implicit real*8 (a-h,o-z)
1478       include 'DIMENSIONS'
1479 #ifdef MPI
1480       include 'mpif.h'
1481 #endif
1482       include 'COMMON.SBRIDGE'
1483       include 'COMMON.IOUNITS'
1484       include 'COMMON.SETUP'
1485       include 'COMMON.CONTROL'
1486       include 'COMMON.MD'
1487       include 'COMMON.INTERACT'
1488 cd      write(iout,*)"homology_partition: lim_odl=",lim_odl,
1489 cd     &   " lim_dih",lim_dih
1490 #ifdef MPI
1491       if (me.eq.king .or. .not. out1file) write (iout,*) "MPI"
1492       call int_bounds(lim_odl,link_start_homo,link_end_homo)
1493       call int_bounds(lim_dih,idihconstr_start_homo,
1494      &  idihconstr_end_homo)
1495       idihconstr_start_homo=idihconstr_start_homo+nnt-1+3
1496       idihconstr_end_homo=idihconstr_end_homo+nnt-1+3
1497       if (me.eq.king .or. .not. out1file)
1498      &  write (iout,*) 'Processor',fg_rank,' CG group',kolor,
1499      &  ' absolute rank',MyRank,
1500      &  ' lim_odl',lim_odl,' link_start=',link_start_homo,
1501      &  ' link_end',link_end_homo,' lim_dih',lim_dih,
1502      &  ' idihconstr_start_homo',idihconstr_start_homo,
1503      &  ' idihconstr_end_homo',idihconstr_end_homo
1504 #else
1505       write (iout,*) "Not MPI"
1506       link_start_homo=1
1507       link_end_homo=lim_odl
1508       idihconstr_start_homo=nnt+3
1509       idihconstr_end_homo=lim_dih+nnt-1+3
1510       write (iout,*)
1511      &  ' lim_odl',lim_odl,' link_start=',link_start_homo,
1512      &  ' link_end',link_end_homo,' lim_dih',lim_dih,
1513      &  ' idihconstr_start_homo',idihconstr_start_homo,
1514      &  ' idihconstr_end_homo',idihconstr_end_homo
1515 #endif
1516       return
1517       end