source/wham/src-M-Pawel/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       enddo
 188 C
 189 C Initialize timing.
 190 C
 191       call set_timers
 192 C
 193 C Initialize variables used in minimization.
 194 C
 195 c     maxfun=5000
 196 c     maxit=2000
 197       maxfun=500
 198       maxit=200
 199       tolf=1.0D-2
 200       rtolf=5.0D-4
 201 C
 202 C Initialize the variables responsible for the mode of gradient storage.
 203 C
 204       nfl=0
 205       icg=1
 206       do i=1,14
 207         do j=1,14
 208           if (print_order(i).eq.j) then
 209             iw(print_order(i))=j
 210             goto 1121
 211           endif
 212         enddo
 213 1121    continue
 214       enddo
 215       calc_grad=.false.
 216 C Set timers and counters for the respective routines
 217       t_func = 0.0d0
 218       t_grad = 0.0d0
 219       t_fhel = 0.0d0
 220       t_fbet = 0.0d0
 221       t_ghel = 0.0d0
 222       t_gbet = 0.0d0
 223       t_viol = 0.0d0
 224       t_gviol = 0.0d0
 225       n_func = 0
 226       n_grad = 0
 227       n_fhel = 0
 228       n_fbet = 0
 229       n_ghel = 0
 230       n_gbet = 0
 231       n_viol = 0
 232       n_gviol = 0
 233       n_map = 0
 234 #ifndef SPLITELE
 235       nprint_ene=nprint_ene-1
 236 #endif
 237       return
 238       end
 239 c-------------------------------------------------------------------------
 240       block data nazwy
 241       implicit real*8 (a-h,o-z)
 242       include 'DIMENSIONS'
 243       include 'DIMENSIONS.ZSCOPT'
 244       include 'COMMON.NAMES'
 245       include 'COMMON.WEIGHTS'
 246       include 'COMMON.FFIELD'
 247       data restyp /
 248      &'DD','DAU','DAI','DDB','DSM','DPR','DLY','DAR','DHI','DAS','DGL',
 249      & 'DSG','DGN','DSN','DTH',
 250      &'DYY','DAL','DTY','DTR','DVA','DLE','DIL','DPN','MED','DCY','ZER',
 251      &'CYS','MET','PHE','ILE','LEU','VAL','TRP','TYR','ALA','GLY','THR',
 252      &'SER','GLN','ASN','GLU','ASP','HIS','ARG','LYS','PRO','SME','DBZ',
 253      &'AIB','ABU','D'/
 254       data onelet /
 255      &'z','z','z','z','z','p','k','r','h','d','e','n','q','s','t','g',
 256      &'a','y','w','v','l','i','f','m','c','x',
 257      &'C','M','F','I','L','V','W','Y','A','G','T',
 258      &'S','Q','N','E','D','H','R','K','P','z','z','z','z','X'/
 259       data potname /'LJ','LJK','BP','GB','GBV'/
 260       data ename /
 261      &   "EVDW SC-SC","EVDW2 SC-p","EES p-p","ECORR4 ","ECORR5 ",
 262      &   "ECORR6 ","EELLO ","ETURN3 ","ETURN4 ","ETURN6 ",
 263      &   "EBE bend","ESC SCloc","ETORS ","ETORSD ","EHPB","EVDWPP",
 264      &   "EVDW2_14","ESTR","ESCCOR","EDIHC","EVDW_T"/
 265       data wname /
 266      &   "WSC","WSCP","WELEC","WCORR","WCORR5","WCORR6","WEL_LOC",
 267      &   "WTURN3","WTURN4","WTURN6","WANG","WSCLOC","WTOR","WTORD",
 268      &   "WHPB","WVDWPP","WSCP14","WBOND","WSCCOR","WDIHC","WSC"/
 269       data ww0 /1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,
 270      &    1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,0.4d0,1.0d0,1.0d0,
 271      &    0.0d0,0.0/
 272       data nprint_ene /21/
 273       data print_order /1,2,3,18,11,12,13,14,4,5,6,7,8,9,10,19,
 274      &  16,15,17,20,21/
 275       end
 276 c---------------------------------------------------------------------------
 277       subroutine init_int_table
 278       implicit real*8 (a-h,o-z)
 279       include 'DIMENSIONS'
 280       include 'DIMENSIONS.ZSCOPT'
 281 #ifdef MPI
 282       include 'mpif.h'
 283 #endif
 284 #ifdef MP
 285       include 'COMMON.INFO'
 286 #endif
 287       include 'COMMON.CHAIN'
 288       include 'COMMON.INTERACT'
 289       include 'COMMON.LOCAL'
 290       include 'COMMON.SBRIDGE'
 291       include 'COMMON.IOUNITS'
 292       logical scheck,lprint
 293 #ifdef MPL
 294       integer my_sc_int(0:max_fg_Procs-1),my_sc_intt(0:max_fg_Procs),
 295      & my_ele_int(0:max_fg_Procs-1),my_ele_intt(0:max_fg_Procs)
 296 C... Determine the numbers of start and end SC-SC interaction
 297 C... to deal with by current processor.
 298       lprint=.false.
 299       if (lprint)
 300      &write (iout,*) 'INIT_INT_TABLE nres=',nres,' nnt=',nnt,' nct=',nct
 301       n_sc_int_tot=(nct-nnt+1)*(nct-nnt)/2-nss
 302       MyRank=MyID-(MyGroup-1)*fgProcs
 303       call int_bounds(n_sc_int_tot,my_sc_inds,my_sc_inde)
 304       if (lprint)
 305      &  write (iout,*) 'Processor',MyID,' MyRank',MyRank,
 306      &  ' n_sc_int_tot',n_sc_int_tot,' my_sc_inds=',my_sc_inds,
 307      &  ' my_sc_inde',my_sc_inde
 308       ind_sctint=0
 309       iatsc_s=0
 310       iatsc_e=0
 311 #endif
 312       lprint=.false.
 313       do i=1,maxres
 314         nint_gr(i)=0
 315         nscp_gr(i)=0
 316         do j=1,maxint_gr
 317           istart(i,1)=0
 318           iend(i,1)=0
 319           ielstart(i)=0
 320           ielend(i)=0
 321           iscpstart(i,1)=0
 322           iscpend(i,1)=0
 323         enddo
 324       enddo
 325       ind_scint=0
 326       ind_scint_old=0
 327 cd    write (iout,*) 'ns=',ns,' nss=',nss,' ihpb,jhpb',
 328 cd   &   (ihpb(i),jhpb(i),i=1,nss)
 329       do i=nnt,nct-1
 330         scheck=.false.
 331         do ii=1,nss
 332           if (ihpb(ii).eq.i+nres) then
 333             scheck=.true.
 334             jj=jhpb(ii)-nres
 335             goto 10
 336           endif
 337         enddo
 338    10   continue
 339 cd      write (iout,*) 'i=',i,' scheck=',scheck,' jj=',jj
 340         if (scheck) then
 341           if (jj.eq.i+1) then
 342 #ifdef MPL
 343             write (iout,*) 'jj=i+1'
 344             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
 345      & iatsc_s,iatsc_e,i+2,nct,nint_gr(i),istart(i,1),iend(i,1),*12)
 346 #else
 347             nint_gr(i)=1
 348             istart(i,1)=i+2
 349             iend(i,1)=nct
 350 #endif
 351           else if (jj.eq.nct) then
 352 #ifdef MPL
 353             write (iout,*) 'jj=nct'
 354             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
 355      &  iatsc_s,iatsc_e,i+1,nct-1,nint_gr(i),istart(i,1),iend(i,1),*12)
 356 #else
 357             nint_gr(i)=1
 358             istart(i,1)=i+1
 359             iend(i,1)=nct-1
 360 #endif
 361           else
 362 #ifdef MPL
 363             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
 364      & iatsc_s,iatsc_e,i+1,jj-1,nint_gr(i),istart(i,1),iend(i,1),*12)
 365             ii=nint_gr(i)+1
 366             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
 367      & iatsc_s,iatsc_e,jj+1,nct,nint_gr(i),istart(i,ii),iend(i,ii),*12)
 368 #else
 369             nint_gr(i)=2
 370             istart(i,1)=i+1
 371             iend(i,1)=jj-1
 372             istart(i,2)=jj+1
 373             iend(i,2)=nct
 374 #endif
 375           endif
 376         else
 377 #ifdef MPL
 378           call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
 379      &    iatsc_s,iatsc_e,i+1,nct,nint_gr(i),istart(i,1),iend(i,1),*12)
 380 #else
 381           nint_gr(i)=1
 382           istart(i,1)=i+1
 383           iend(i,1)=nct
 384           ind_scint=int_scint+nct-i
 385 #endif
 386         endif
 387 #ifdef MPL
 388         ind_scint_old=ind_scint
 389 #endif
 390       enddo
 391    12 continue
 392 #ifndef MPL
 393       iatsc_s=nnt
 394       iatsc_e=nct-1
 395 #endif
 396 #ifdef MPL
 397       if (lprint) then
 398         write (iout,*) 'Processor',MyID,' Group',MyGroup
 399         write (iout,*) 'iatsc_s=',iatsc_s,' iatsc_e=',iatsc_e
 400       endif
 401 #endif
 402       if (lprint) then
 403       write (iout,'(a)') 'Interaction array:'
 404       do i=iatsc_s,iatsc_e
 405         write (iout,'(i3,2(2x,2i3))')
 406      & i,(istart(i,iint),iend(i,iint),iint=1,nint_gr(i))
 407       enddo
 408       endif
 409       ispp=2
 410 #ifdef MPL
 411 C Now partition the electrostatic-interaction array
 412       npept=nct-nnt
 413       nele_int_tot=(npept-ispp)*(npept-ispp+1)/2
 414       call int_bounds(nele_int_tot,my_ele_inds,my_ele_inde)
 415       if (lprint)
 416      & write (iout,*) 'Processor',MyID,' MyRank',MyRank,
 417      &  ' nele_int_tot',nele_int_tot,' my_ele_inds=',my_ele_inds,
 418      &               ' my_ele_inde',my_ele_inde
 419       iatel_s=0
 420       iatel_e=0
 421       ind_eleint=0
 422       ind_eleint_old=0
 423       do i=nnt,nct-3
 424         ijunk=0
 425         call int_partition(ind_eleint,my_ele_inds,my_ele_inde,i,
 426      &    iatel_s,iatel_e,i+ispp,nct-1,ijunk,ielstart(i),ielend(i),*13)
 427       enddo ! i
 428    13 continue
 429 #else
 430       iatel_s=nnt
 431       iatel_e=nct-3
 432       do i=iatel_s,iatel_e
 433         ielstart(i)=i+2
 434         ielend(i)=nct-1
 435       enddo
 436 #endif
 437       if (lprint) then
 438         write (iout,'(a)') 'Electrostatic interaction array:'
 439         do i=iatel_s,iatel_e
 440           write (iout,'(i3,2(2x,2i3))') i,ielstart(i),ielend(i)
 441         enddo
 442       endif ! lprint
 443 c     iscp=3
 444       iscp=2
 445 C Partition the SC-p interaction array
 446 #ifdef MPL
 447       nscp_int_tot=(npept-iscp+1)*(npept-iscp+1)
 448       call int_bounds(nscp_int_tot,my_scp_inds,my_scp_inde)
 449       if (lprint)
 450      & write (iout,*) 'Processor',MyID,' MyRank',MyRank,
 451      &  ' nscp_int_tot',nscp_int_tot,' my_scp_inds=',my_scp_inds,
 452      &               ' my_scp_inde',my_scp_inde
 453       iatscp_s=0
 454       iatscp_e=0
 455       ind_scpint=0
 456       ind_scpint_old=0
 457       do i=nnt,nct-1
 458         if (i.lt.nnt+iscp) then
 459 cd        write (iout,*) 'i.le.nnt+iscp'
 460           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
 461      &      iatscp_s,iatscp_e,i+iscp,nct,nscp_gr(i),iscpstart(i,1),
 462      &      iscpend(i,1),*14)
 463         else if (i.gt.nct-iscp) then
 464 cd        write (iout,*) 'i.gt.nct-iscp'
 465           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
 466      &      iatscp_s,iatscp_e,nnt,i-iscp,nscp_gr(i),iscpstart(i,1),
 467      &      iscpend(i,1),*14)
 468         else
 469           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
 470      &      iatscp_s,iatscp_e,nnt,i-iscp,nscp_gr(i),iscpstart(i,1),
 471      &      iscpend(i,1),*14)
 472           ii=nscp_gr(i)+1
 473           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
 474      &      iatscp_s,iatscp_e,i+iscp,nct,nscp_gr(i),iscpstart(i,ii),
 475      &      iscpend(i,ii),*14)
 476         endif
 477       enddo ! i
 478    14 continue
 479 #else
 480       iatscp_s=nnt
 481       iatscp_e=nct-1
 482       do i=nnt,nct-1
 483         if (i.lt.nnt+iscp) then
 484           nscp_gr(i)=1
 485           iscpstart(i,1)=i+iscp
 486           iscpend(i,1)=nct
 487         elseif (i.gt.nct-iscp) then
 488           nscp_gr(i)=1
 489           iscpstart(i,1)=nnt
 490           iscpend(i,1)=i-iscp
 491         else
 492           nscp_gr(i)=2
 493           iscpstart(i,1)=nnt
 494           iscpend(i,1)=i-iscp
 495           iscpstart(i,2)=i+iscp
 496           iscpend(i,2)=nct
 497         endif
 498       enddo ! i
 499 #endif
 500       if (lprint) then
 501         write (iout,'(a)') 'SC-p interaction array:'
 502         do i=iatscp_s,iatscp_e
 503           write (iout,'(i3,2(2x,2i3))')
 504      &         i,(iscpstart(i,j),iscpend(i,j),j=1,nscp_gr(i))
 505         enddo
 506       endif ! lprint
 507 C Partition local interactions
 508 #ifdef MPL
 509       call int_bounds(nres-2,loc_start,loc_end)
 510       loc_start=loc_start+1
 511       loc_end=loc_end+1
 512       call int_bounds(nres-2,ithet_start,ithet_end)
 513       ithet_start=ithet_start+2
 514       ithet_end=ithet_end+2
 515       call int_bounds(nct-nnt-2,iphi_start,iphi_end)
 516       iphi_start=iphi_start+nnt+2
 517       iphi_end=iphi_end+nnt+2
 518       call int_bounds(nres-3,itau_start,itau_end)
 519       itau_start=itau_start+3
 520       itau_end=itau_end+3
 521       if (lprint) then
 522         write (iout,*) 'Processor:',MyID,
 523      & ' loc_start',loc_start,' loc_end',loc_end,
 524      & ' ithet_start',ithet_start,' ithet_end',ithet_end,
 525      & ' iphi_start',iphi_start,' iphi_end',iphi_end
 526         write (*,*) 'Processor:',MyID,
 527      & ' loc_start',loc_start,' loc_end',loc_end,
 528      & ' ithet_start',ithet_start,' ithet_end',ithet_end,
 529      & ' iphi_start',iphi_start,' iphi_end',iphi_end
 530       endif
 531       if (fgprocs.gt.1 .and. MyID.eq.BossID) then
 532         write(iout,'(i10,a,i10,a,i10,a/a,i3,a)') n_sc_int_tot,' SC-SC ',
 533      & nele_int_tot,' electrostatic and ',nscp_int_tot,
 534      & ' SC-p interactions','were distributed among',fgprocs,
 535      & ' fine-grain processors.'
 536       endif
 537 #else
 538       loc_start=2
 539       loc_end=nres-1
 540       ithet_start=3
 541       ithet_end=nres
 542       iphi_start=nnt+3
 543       iphi_end=nct
 544       itau_start=4
 545       itau_end=nres
 546 #endif
 547       return
 548       end
 549 c---------------------------------------------------------------------------
 550       subroutine int_partition(int_index,lower_index,upper_index,atom,
 551      & at_start,at_end,first_atom,last_atom,int_gr,jat_start,jat_end,*)
 552       implicit real*8 (a-h,o-z)
 553       include 'DIMENSIONS'
 554       include 'COMMON.IOUNITS'
 555       integer int_index,lower_index,upper_index,atom,at_start,at_end,
 556      & first_atom,last_atom,int_gr,jat_start,jat_end
 557       logical lprn
 558       lprn=.false.
 559       if (lprn) write (iout,*) 'int_index=',int_index
 560       int_index_old=int_index
 561       int_index=int_index+last_atom-first_atom+1
 562       if (lprn)
 563      &   write (iout,*) 'int_index=',int_index,
 564      &               ' int_index_old',int_index_old,
 565      &               ' lower_index=',lower_index,
 566      &               ' upper_index=',upper_index,
 567      &               ' atom=',atom,' first_atom=',first_atom,
 568      &               ' last_atom=',last_atom
 569       if (int_index.ge.lower_index) then
 570         int_gr=int_gr+1
 571         if (at_start.eq.0) then
 572           at_start=atom
 573           jat_start=first_atom-1+lower_index-int_index_old
 574         else
 575           jat_start=first_atom
 576         endif
 577         if (lprn) write (iout,*) 'jat_start',jat_start
 578         if (int_index.ge.upper_index) then
 579           at_end=atom
 580           jat_end=first_atom-1+upper_index-int_index_old
 581           return1
 582         else
 583           jat_end=last_atom
 584         endif
 585         if (lprn) write (iout,*) 'jat_end',jat_end
 586       endif
 587       return
 588       end
 589 c------------------------------------------------------------------------------
 590       subroutine hpb_partition
 591       implicit real*8 (a-h,o-z)
 592       include 'DIMENSIONS'
 593       include 'COMMON.SBRIDGE'
 594       include 'COMMON.IOUNITS'
 595 #ifdef MPL
 596       include 'COMMON.INFO'
 597       call int_bounds(nhpb,link_start,link_end)
 598 #else
 599       link_start=1
 600       link_end=nhpb
 601 #endif
 602 cd    write (iout,*) 'Processor',MyID,' MyRank',MyRank,
 603 cd   &  ' nhpb',nhpb,' link_start=',link_start,
 604 cd   &  ' link_end',link_end
 605       return
 606       end