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