source/cluster/wham/src-M/initialize_p.F

   1       subroutine initialize
   2 C
   3 C Define constants and zero out tables.
   4 C
   5       implicit real*8 (a-h,o-z)
   6       include 'DIMENSIONS'
   7       include 'sizesclu.dat'
   8 #ifdef MPI
   9       include 'mpif.h'
  10 #endif
  11       include 'COMMON.IOUNITS'
  12       include 'COMMON.CHAIN'
  13       include 'COMMON.INTERACT'
  14       include 'COMMON.GEO'
  15       include 'COMMON.LOCAL'
  16       include 'COMMON.TORSION'
  17       include 'COMMON.FFIELD'
  18       include 'COMMON.SBRIDGE'
  19       include 'COMMON.MINIM'
  20       include 'COMMON.DERIV'
  21       include "COMMON.NAMES"
  22       include "COMMON.TIME1"
  23 C
  24 C The following is just to define auxiliary variables used in angle conversion
  25 C
  26       pi=4.0D0*datan(1.0D0)
  27       dwapi=2.0D0*pi
  28       dwapi3=dwapi/3.0D0
  29       pipol=0.5D0*pi
  30       deg2rad=pi/180.0D0
  31       rad2deg=1.0D0/deg2rad
  32       angmin=10.0D0*deg2rad
  33       Rgas = 1.987D-3
  34 C
  35 C Define I/O units.
  36 C
  37       inp=    1
  38       iout=   2
  39       ipdbin= 3
  40       ipdb=   7
  41       imol2= 18
  42       jplot= 19
  43       jstatin=10
  44       imol2=  4
  45       igeom=  8
  46       intin=  9
  47       ithep= 11
  48       irotam=12
  49       itorp= 13
  50       itordp= 23
  51       ielep= 14
  52       isidep=15
  53       isidep1=22
  54       iscpp=25
  55       icbase=16
  56       ifourier=20
  57       istat= 17
  58       ibond=28
  59       isccor=29
  60       jrms=30
  61 C
  62 C Set default weights of the energy terms.
  63 C
  64       wlong=1.0D0
  65       welec=1.0D0
  66       wtor =1.0D0
  67       wang =1.0D0
  68       wscloc=1.0D0
  69       wstrain=1.0D0
  70 C
  71 C Zero out tables.
  72 C
  73       ndih_constr=0
  74       do i=1,maxres2
  75         do j=1,3
  76           c(j,i)=0.0D0
  77           dc(j,i)=0.0D0
  78         enddo
  79       enddo
  80       do i=1,maxres
  81         do j=1,3
  82           xloc(j,i)=0.0D0
  83         enddo
  84       enddo
  85       do i=1,ntyp
  86         do j=1,ntyp
  87           aa(i,j)=0.0D0
  88           bb(i,j)=0.0D0
  89           augm(i,j)=0.0D0
  90           sigma(i,j)=0.0D0
  91           r0(i,j)=0.0D0
  92           chi(i,j)=0.0D0
  93         enddo
  94         do j=1,2
  95           bad(i,j)=0.0D0
  96         enddo
  97         chip(i)=0.0D0
  98         alp(i)=0.0D0
  99         sigma0(i)=0.0D0
 100         sigii(i)=0.0D0
 101         rr0(i)=0.0D0
 102         a0thet(i)=0.0D0
 103         do j=1,2
 104          do ichir1=-1,1
 105           do ichir2=-1,1
 106           athet(j,i,ichir1,ichir2)=0.0D0
 107           bthet(j,i,ichir1,ichir2)=0.0D0
 108           enddo
 109          enddo
 110         enddo
 111         do j=0,3
 112           polthet(j,i)=0.0D0
 113         enddo
 114         do j=1,3
 115           gthet(j,i)=0.0D0
 116         enddo
 117         theta0(i)=0.0D0
 118         sig0(i)=0.0D0
 119         sigc0(i)=0.0D0
 120         do j=1,maxlob
 121           bsc(j,i)=0.0D0
 122           do k=1,3
 123             censc(k,j,i)=0.0D0
 124           enddo
 125           do k=1,3
 126             do l=1,3
 127               gaussc(l,k,j,i)=0.0D0
 128             enddo
 129           enddo
 130           nlob(i)=0
 131         enddo
 132       enddo
 133       nlob(ntyp1)=0
 134       dsc(ntyp1)=0.0D0
 135       do i=-maxtor,maxtor
 136         itortyp(i)=0
 137        do iblock=1,2
 138         do j=-maxtor,maxtor
 139           do k=1,maxterm
 140             v1(k,j,i,iblock)=0.0D0
 141             v2(k,j,i,iblock)=0.0D0
 142            enddo
 143          enddo
 144        enddo
 145       enddo
 146       do iblock=1,2
 147        do i=-maxtor,maxtor
 148         do j=-maxtor,maxtor
 149          do k=-maxtor,maxtor
 150           do l=1,maxtermd_1
 151             v1c(1,l,i,j,k,iblock)=0.0D0
 152             v1s(1,l,i,j,k,iblock)=0.0D0
 153             v1c(2,l,i,j,k,iblock)=0.0D0
 154             v1s(2,l,i,j,k,iblock)=0.0D0
 155           enddo !l
 156           do l=1,maxtermd_2
 157            do m=1,maxtermd_2
 158             v2c(m,l,i,j,k,iblock)=0.0D0
 159             v2s(m,l,i,j,k,iblock)=0.0D0
 160            enddo !m
 161           enddo !l
 162         enddo !k
 163        enddo !j
 164       enddo !i
 165       enddo !iblock
 166       do i=1,maxres
 167         itype(i)=0
 168         itel(i)=0
 169       enddo
 170 C Initialize the bridge arrays
 171       ns=0
 172       nss=0
 173       nhpb=0
 174       do i=1,maxss
 175         iss(i)=0
 176       enddo
 177       do i=1,maxss
 178         dhpb(i)=0.0D0
 179       enddo
 180       do i=1,maxss
 181         ihpb(i)=0
 182         jhpb(i)=0
 183       enddo
 184 C
 185 C Initialize timing.
 186 C
 187       call set_timers
 188 C
 189 C Initialize variables used in minimization.
 190 C
 191 c     maxfun=5000
 192 c     maxit=2000
 193       maxfun=500
 194       maxit=200
 195       tolf=1.0D-2
 196       rtolf=5.0D-4
 197 C
 198 C Initialize the variables responsible for the mode of gradient storage.
 199 C
 200       nfl=0
 201       icg=1
 202       do i=1,14
 203         do j=1,14
 204           if (print_order(i).eq.j) then
 205             iw(print_order(i))=j
 206             goto 1121
 207           endif
 208         enddo
 209 1121    continue
 210       enddo
 211       calc_grad=.false.
 212 C Set timers and counters for the respective routines
 213       t_func = 0.0d0
 214       t_grad = 0.0d0
 215       t_fhel = 0.0d0
 216       t_fbet = 0.0d0
 217       t_ghel = 0.0d0
 218       t_gbet = 0.0d0
 219       t_viol = 0.0d0
 220       t_gviol = 0.0d0
 221       n_func = 0
 222       n_grad = 0
 223       n_fhel = 0
 224       n_fbet = 0
 225       n_ghel = 0
 226       n_gbet = 0
 227       n_viol = 0
 228       n_gviol = 0
 229       n_map = 0
 230 #ifndef SPLITELE
 231       nprint_ene=nprint_ene-1
 232 #endif
 233       return
 234       end
 235 c-------------------------------------------------------------------------
 236       block data nazwy
 237       implicit real*8 (a-h,o-z)
 238       include 'DIMENSIONS'
 239       include 'sizesclu.dat'
 240       include 'COMMON.NAMES'
 241       include 'COMMON.FFIELD'
 242       data restyp /
 243      &'DD','DAU','DAI','DDB','DSM','DPR','DLY','DAR','DHI','DAS','DGL',
 244      & 'DSG','DGN','DSN','DTH',
 245      &'DYY','DAL','DTY','DTR','DVA','DLE','DIL','DPN','MED','DCY','ZER',
 246      &'CYS','MET','PHE','ILE','LEU','VAL','TRP','TYR','ALA','GLY','THR',
 247      &'SER','GLN','ASN','GLU','ASP','HIS','ARG','LYS','PRO','SME','DBZ',
 248      &'AIB','ABU','D'/
 249       data onelet /
 250      &'z','z','z','z','z','p','k','r','h','d','e','n','q','s','t','g',
 251      &'a','y','w','v','l','i','f','m','c','x',
 252      &'C','M','F','I','L','V','W','Y','A','G','T',
 253      &'S','Q','N','E','D','H','R','K','P','z','z','z','z','X'/
 254       data potname /'LJ','LJK','BP','GB','GBV'/
 255       data ename /
 256      &   "EVDW SC-SC","EVDW2 SC-p","EES p-p","ECORR4 ","ECORR5 ",
 257      &   "ECORR6 ","EELLO ","ETURN3 ","ETURN4 ","ETURN6 ",
 258      &   "EBE bend","ESC SCloc","ETORS ","ETORSD ","EHPB","EVDWPP",
 259      &   "ESTR","ESCCOR","EVDW2_14","","EVDW_T"/
 260       data wname /
 261      &   "WSC","WSCP","WELEC","WCORR","WCORR5","WCORR6","WEL_LOC",
 262      &   "WTURN3","WTURN4","WTURN6","WANG","WSCLOC","WTOR","WTORD",
 263      &   "WHPB","WVDWPP","WBOND","WSCCOR","WSCP14","","WSC"/
 264       data nprint_ene /21/
 265       data print_order /1,2,3,18,11,12,13,14,4,5,6,7,8,9,10,19,
 266      &  16,15,17,20,21/
 267       end
 268 c---------------------------------------------------------------------------
 269       subroutine init_int_table
 270       implicit real*8 (a-h,o-z)
 271       include 'DIMENSIONS'
 272       include 'sizesclu.dat'
 273 #ifdef MPI
 274       include 'mpif.h'
 275 #endif
 276 #ifdef MPL
 277       include 'COMMON.INFO'
 278 #endif
 279       include 'COMMON.CHAIN'
 280       include 'COMMON.INTERACT'
 281       include 'COMMON.LOCAL'
 282       include 'COMMON.SBRIDGE'
 283       include 'COMMON.IOUNITS'
 284       logical scheck,lprint
 285 #ifdef MPL
 286       integer my_sc_int(0:max_fg_Procs-1),my_sc_intt(0:max_fg_Procs),
 287      & my_ele_int(0:max_fg_Procs-1),my_ele_intt(0:max_fg_Procs)
 288 C... Determine the numbers of start and end SC-SC interaction
 289 C... to deal with by current processor.
 290       lprint=.false.
 291       if (lprint)
 292      &write (iout,*) 'INIT_INT_TABLE nres=',nres,' nnt=',nnt,' nct=',nct
 293       n_sc_int_tot=(nct-nnt+1)*(nct-nnt)/2-nss
 294       MyRank=MyID-(MyGroup-1)*fgProcs
 295       call int_bounds(n_sc_int_tot,my_sc_inds,my_sc_inde)
 296       if (lprint)
 297      &  write (iout,*) 'Processor',MyID,' MyRank',MyRank,
 298      &  ' n_sc_int_tot',n_sc_int_tot,' my_sc_inds=',my_sc_inds,
 299      &  ' my_sc_inde',my_sc_inde
 300       ind_sctint=0
 301       iatsc_s=0
 302       iatsc_e=0
 303 #endif
 304       lprint=.false.
 305       do i=1,maxres
 306         nint_gr(i)=0
 307         nscp_gr(i)=0
 308         do j=1,maxint_gr
 309           istart(i,1)=0
 310           iend(i,1)=0
 311           ielstart(i)=0
 312           ielend(i)=0
 313           iscpstart(i,1)=0
 314           iscpend(i,1)=0
 315         enddo
 316       enddo
 317       ind_scint=0
 318       ind_scint_old=0
 319 cd    write (iout,*) 'ns=',ns,' nss=',nss,' ihpb,jhpb',
 320 cd   &   (ihpb(i),jhpb(i),i=1,nss)
 321       do i=nnt,nct-1
 322         scheck=.false.
 323         do ii=1,nss
 324           if (ihpb(ii).eq.i+nres) then
 325             scheck=.true.
 326             jj=jhpb(ii)-nres
 327             goto 10
 328           endif
 329         enddo
 330    10   continue
 331 cd      write (iout,*) 'i=',i,' scheck=',scheck,' jj=',jj
 332         if (scheck) then
 333           if (jj.eq.i+1) then
 334 #ifdef MPL
 335             write (iout,*) 'jj=i+1'
 336             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
 337      & iatsc_s,iatsc_e,i+2,nct,nint_gr(i),istart(i,1),iend(i,1),*12)
 338 #else
 339             nint_gr(i)=1
 340             istart(i,1)=i+2
 341             iend(i,1)=nct
 342 #endif
 343           else if (jj.eq.nct) then
 344 #ifdef MPL
 345             write (iout,*) 'jj=nct'
 346             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
 347      &  iatsc_s,iatsc_e,i+1,nct-1,nint_gr(i),istart(i,1),iend(i,1),*12)
 348 #else
 349             nint_gr(i)=1
 350             istart(i,1)=i+1
 351             iend(i,1)=nct-1
 352 #endif
 353           else
 354 #ifdef MPL
 355             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
 356      & iatsc_s,iatsc_e,i+1,jj-1,nint_gr(i),istart(i,1),iend(i,1),*12)
 357             ii=nint_gr(i)+1
 358             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
 359      & iatsc_s,iatsc_e,jj+1,nct,nint_gr(i),istart(i,ii),iend(i,ii),*12)
 360 #else
 361             nint_gr(i)=2
 362             istart(i,1)=i+1
 363             iend(i,1)=jj-1
 364             istart(i,2)=jj+1
 365             iend(i,2)=nct
 366 #endif
 367           endif
 368         else
 369 #ifdef MPL
 370           call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
 371      &    iatsc_s,iatsc_e,i+1,nct,nint_gr(i),istart(i,1),iend(i,1),*12)
 372 #else
 373           nint_gr(i)=1
 374           istart(i,1)=i+1
 375           iend(i,1)=nct
 376           ind_scint=int_scint+nct-i
 377 #endif
 378         endif
 379 #ifdef MPL
 380         ind_scint_old=ind_scint
 381 #endif
 382       enddo
 383    12 continue
 384 #ifndef MPL
 385       iatsc_s=nnt
 386       iatsc_e=nct-1
 387 #endif
 388 #ifdef MPL
 389       if (lprint) then
 390         write (iout,*) 'Processor',MyID,' Group',MyGroup
 391         write (iout,*) 'iatsc_s=',iatsc_s,' iatsc_e=',iatsc_e
 392       endif
 393 #endif
 394       if (lprint) then
 395       write (iout,'(a)') 'Interaction array:'
 396       do i=iatsc_s,iatsc_e
 397         write (iout,'(i3,2(2x,2i3))')
 398      & i,(istart(i,iint),iend(i,iint),iint=1,nint_gr(i))
 399       enddo
 400       endif
 401       ispp=2
 402 #ifdef MPL
 403 C Now partition the electrostatic-interaction array
 404       npept=nct-nnt
 405       nele_int_tot=(npept-ispp)*(npept-ispp+1)/2
 406       call int_bounds(nele_int_tot,my_ele_inds,my_ele_inde)
 407       if (lprint)
 408      & write (iout,*) 'Processor',MyID,' MyRank',MyRank,
 409      &  ' nele_int_tot',nele_int_tot,' my_ele_inds=',my_ele_inds,
 410      &               ' my_ele_inde',my_ele_inde
 411       iatel_s=0
 412       iatel_e=0
 413       ind_eleint=0
 414       ind_eleint_old=0
 415       do i=nnt,nct-3
 416         ijunk=0
 417         call int_partition(ind_eleint,my_ele_inds,my_ele_inde,i,
 418      &    iatel_s,iatel_e,i+ispp,nct-1,ijunk,ielstart(i),ielend(i),*13)
 419       enddo ! i
 420    13 continue
 421 #else
 422       iatel_s=nnt
 423       iatel_e=nct-3
 424       do i=iatel_s,iatel_e
 425         ielstart(i)=i+2
 426         ielend(i)=nct-1
 427       enddo
 428 #endif
 429       if (lprint) then
 430         write (iout,'(a)') 'Electrostatic interaction array:'
 431         do i=iatel_s,iatel_e
 432           write (iout,'(i3,2(2x,2i3))') i,ielstart(i),ielend(i)
 433         enddo
 434       endif ! lprint
 435 c     iscp=3
 436       iscp=2
 437 C Partition the SC-p interaction array
 438 #ifdef MPL
 439       nscp_int_tot=(npept-iscp+1)*(npept-iscp+1)
 440       call int_bounds(nscp_int_tot,my_scp_inds,my_scp_inde)
 441       if (lprint)
 442      & write (iout,*) 'Processor',MyID,' MyRank',MyRank,
 443      &  ' nscp_int_tot',nscp_int_tot,' my_scp_inds=',my_scp_inds,
 444      &               ' my_scp_inde',my_scp_inde
 445       iatscp_s=0
 446       iatscp_e=0
 447       ind_scpint=0
 448       ind_scpint_old=0
 449       do i=nnt,nct-1
 450         if (i.lt.nnt+iscp) then
 451 cd        write (iout,*) 'i.le.nnt+iscp'
 452           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
 453      &      iatscp_s,iatscp_e,i+iscp,nct,nscp_gr(i),iscpstart(i,1),
 454      &      iscpend(i,1),*14)
 455         else if (i.gt.nct-iscp) then
 456 cd        write (iout,*) 'i.gt.nct-iscp'
 457           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
 458      &      iatscp_s,iatscp_e,nnt,i-iscp,nscp_gr(i),iscpstart(i,1),
 459      &      iscpend(i,1),*14)
 460         else
 461           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
 462      &      iatscp_s,iatscp_e,nnt,i-iscp,nscp_gr(i),iscpstart(i,1),
 463      &      iscpend(i,1),*14)
 464           ii=nscp_gr(i)+1
 465           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
 466      &      iatscp_s,iatscp_e,i+iscp,nct,nscp_gr(i),iscpstart(i,ii),
 467      &      iscpend(i,ii),*14)
 468         endif
 469       enddo ! i
 470    14 continue
 471 #else
 472       iatscp_s=nnt
 473       iatscp_e=nct-1
 474       do i=nnt,nct-1
 475         if (i.lt.nnt+iscp) then
 476           nscp_gr(i)=1
 477           iscpstart(i,1)=i+iscp
 478           iscpend(i,1)=nct
 479         elseif (i.gt.nct-iscp) then
 480           nscp_gr(i)=1
 481           iscpstart(i,1)=nnt
 482           iscpend(i,1)=i-iscp
 483         else
 484           nscp_gr(i)=2
 485           iscpstart(i,1)=nnt
 486           iscpend(i,1)=i-iscp
 487           iscpstart(i,2)=i+iscp
 488           iscpend(i,2)=nct
 489         endif
 490       enddo ! i
 491 #endif
 492       if (lprint) then
 493         write (iout,'(a)') 'SC-p interaction array:'
 494         do i=iatscp_s,iatscp_e
 495           write (iout,'(i3,2(2x,2i3))')
 496      &         i,(iscpstart(i,j),iscpend(i,j),j=1,nscp_gr(i))
 497         enddo
 498       endif ! lprint
 499 C Partition local interactions
 500 #ifdef MPL
 501       call int_bounds(nres-2,loc_start,loc_end)
 502       loc_start=loc_start+1
 503       loc_end=loc_end+1
 504       call int_bounds(nres-2,ithet_start,ithet_end)
 505       ithet_start=ithet_start+2
 506       ithet_end=ithet_end+2
 507       call int_bounds(nct-nnt-2,iphi_start,iphi_end)
 508       iphi_start=iphi_start+nnt+2
 509       iphi_end=iphi_end+nnt+2
 510       call int_bounds(nres-3,itau_start,itau_end)
 511       itau_start=itau_start+3
 512       itau_end=itau_end+3
 513       if (lprint) then
 514         write (iout,*) 'Processor:',MyID,
 515      & ' loc_start',loc_start,' loc_end',loc_end,
 516      & ' ithet_start',ithet_start,' ithet_end',ithet_end,
 517      & ' iphi_start',iphi_start,' iphi_end',iphi_end
 518         write (*,*) 'Processor:',MyID,
 519      & ' loc_start',loc_start,' loc_end',loc_end,
 520      & ' ithet_start',ithet_start,' ithet_end',ithet_end,
 521      & ' iphi_start',iphi_start,' iphi_end',iphi_end
 522       endif
 523       if (fgprocs.gt.1 .and. MyID.eq.BossID) then
 524         write(iout,'(i10,a,i10,a,i10,a/a,i3,a)') n_sc_int_tot,' SC-SC ',
 525      & nele_int_tot,' electrostatic and ',nscp_int_tot,
 526      & ' SC-p interactions','were distributed among',fgprocs,
 527      & ' fine-grain processors.'
 528       endif
 529 #else
 530       loc_start=2
 531       loc_end=nres-1
 532       ithet_start=3
 533       ithet_end=nres
 534       iphi_start=nnt+3
 535       iphi_end=nct
 536       itau_start=4
 537       itau_end=nres
 538 #endif
 539       return
 540       end
 541 c---------------------------------------------------------------------------
 542       subroutine int_partition(int_index,lower_index,upper_index,atom,
 543      & at_start,at_end,first_atom,last_atom,int_gr,jat_start,jat_end,*)
 544       implicit real*8 (a-h,o-z)
 545       include 'DIMENSIONS'
 546       include 'COMMON.IOUNITS'
 547       integer int_index,lower_index,upper_index,atom,at_start,at_end,
 548      & first_atom,last_atom,int_gr,jat_start,jat_end
 549       logical lprn
 550       lprn=.false.
 551       if (lprn) write (iout,*) 'int_index=',int_index
 552       int_index_old=int_index
 553       int_index=int_index+last_atom-first_atom+1
 554       if (lprn)
 555      &   write (iout,*) 'int_index=',int_index,
 556      &               ' int_index_old',int_index_old,
 557      &               ' lower_index=',lower_index,
 558      &               ' upper_index=',upper_index,
 559      &               ' atom=',atom,' first_atom=',first_atom,
 560      &               ' last_atom=',last_atom
 561       if (int_index.ge.lower_index) then
 562         int_gr=int_gr+1
 563         if (at_start.eq.0) then
 564           at_start=atom
 565           jat_start=first_atom-1+lower_index-int_index_old
 566         else
 567           jat_start=first_atom
 568         endif
 569         if (lprn) write (iout,*) 'jat_start',jat_start
 570         if (int_index.ge.upper_index) then
 571           at_end=atom
 572           jat_end=first_atom-1+upper_index-int_index_old
 573           return1
 574         else
 575           jat_end=last_atom
 576         endif
 577         if (lprn) write (iout,*) 'jat_end',jat_end
 578       endif
 579       return
 580       end
 581 c------------------------------------------------------------------------------
 582       subroutine hpb_partition
 583       implicit real*8 (a-h,o-z)
 584       include 'DIMENSIONS'
 585       include 'COMMON.SBRIDGE'
 586       include 'COMMON.IOUNITS'
 587 #ifdef MPL
 588       include 'COMMON.INFO'
 589       call int_bounds(nhpb,link_start,link_end)
 590 #else
 591       link_start=1
 592       link_end=nhpb
 593 #endif
 594 cd    write (iout,*) 'Processor',MyID,' MyRank',MyRank,
 595 cd   &  ' nhpb',nhpb,' link_start=',link_start,
 596 cd   &  ' link_end',link_end
 597       return
 598       end