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