source/wham/src/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       ithep_pdb=51
  46       irotam=12
  47       irotam_pdb=52
  48       itorp= 13
  49       itordp= 23
  50       ielep= 14
  51       isidep=15
  52       isidep1=22
  53       iscpp=25
  54       icbase=16
  55       ifourier=20
  56       istat= 17
  57       ientin=18
  58       ientout=19
  59       ibond=28
  60       isccor=29
  61 C
  62 C WHAM files
  63 C
  64       ihist=30
  65       iweight=31
  66       izsc=32
  67 C
  68 C Set default weights of the energy terms.
  69 C
  70       wlong=1.0D0
  71       welec=1.0D0
  72       wtor =1.0D0
  73       wang =1.0D0
  74       wscloc=1.0D0
  75       wstrain=1.0D0
  76 C
  77 C Zero out tables.
  78 C
  79       ndih_constr=0
  80       do i=1,maxres2
  81         do j=1,3
  82           c(j,i)=0.0D0
  83           dc(j,i)=0.0D0
  84         enddo
  85       enddo
  86       do i=1,maxres
  87         do j=1,3
  88           xloc(j,i)=0.0D0
  89         enddo
  90       enddo
  91       do i=1,ntyp
  92         do j=1,ntyp
  93           aa(i,j)=0.0D0
  94           bb(i,j)=0.0D0
  95           augm(i,j)=0.0D0
  96           sigma(i,j)=0.0D0
  97           r0(i,j)=0.0D0
  98           chi(i,j)=0.0D0
  99         enddo
 100         do j=1,2
 101           bad(i,j)=0.0D0
 102         enddo
 103         chip(i)=0.0D0
 104         alp(i)=0.0D0
 105         sigma0(i)=0.0D0
 106         sigii(i)=0.0D0
 107         rr0(i)=0.0D0
 108         a0thet(i)=0.0D0
 109         do j=1,2
 110          do ichir1=-1,1
 111           do ichir2=-1,1
 112           athet(j,i,ichir1,ichir2)=0.0D0
 113           bthet(j,i,ichir1,ichir2)=0.0D0
 114           enddo
 115          enddo
 116         enddo
 117         do j=0,3
 118           polthet(j,i)=0.0D0
 119         enddo
 120         do j=1,3
 121           gthet(j,i)=0.0D0
 122         enddo
 123         theta0(i)=0.0D0
 124         sig0(i)=0.0D0
 125         sigc0(i)=0.0D0
 126         do j=1,maxlob
 127           bsc(j,i)=0.0D0
 128           do k=1,3
 129             censc(k,j,i)=0.0D0
 130           enddo
 131           do k=1,3
 132             do l=1,3
 133               gaussc(l,k,j,i)=0.0D0
 134             enddo
 135           enddo
 136           nlob(i)=0
 137         enddo
 138       enddo
 139       nlob(ntyp1)=0
 140       dsc(ntyp1)=0.0D0
 141       do i=-maxtor,maxtor
 142         itortyp(i)=0
 143        do iblock=1,2
 144         do j=-maxtor,maxtor
 145           do k=1,maxterm
 146             v1(k,j,i,iblock)=0.0D0
 147             v2(k,j,i,iblock)=0.0D0
 148           enddo
 149         enddo
 150         enddo
 151       enddo
 152       do iblock=1,2
 153        do i=-maxtor,maxtor
 154         do j=-maxtor,maxtor
 155          do k=-maxtor,maxtor
 156           do l=1,maxtermd_1
 157             v1c(1,l,i,j,k,iblock)=0.0D0
 158             v1s(1,l,i,j,k,iblock)=0.0D0
 159             v1c(2,l,i,j,k,iblock)=0.0D0
 160             v1s(2,l,i,j,k,iblock)=0.0D0
 161           enddo !l
 162           do l=1,maxtermd_2
 163            do m=1,maxtermd_2
 164             v2c(m,l,i,j,k,iblock)=0.0D0
 165             v2s(m,l,i,j,k,iblock)=0.0D0
 166            enddo !m
 167           enddo !l
 168         enddo !k
 169        enddo !j
 170       enddo !i
 171       enddo !iblock
 172       do i=1,maxres
 173         itype(i)=0
 174         itel(i)=0
 175       enddo
 176 C Initialize the bridge arrays
 177       ns=0
 178       nss=0
 179       nhpb=0
 180       do i=1,maxss
 181         iss(i)=0
 182       enddo
 183       do i=1,maxdim
 184         dhpb(i)=0.0D0
 185       enddo
 186       do i=1,maxres
 187         ihpb(i)=0
 188         jhpb(i)=0
 189       enddo
 190 C
 191 C Initialize timing.
 192 C
 193       call set_timers
 194 C
 195 C Initialize variables used in minimization.
 196 C
 197 c     maxfun=5000
 198 c     maxit=2000
 199       maxfun=500
 200       maxit=200
 201       tolf=1.0D-2
 202       rtolf=5.0D-4
 203 C
 204 C Initialize the variables responsible for the mode of gradient storage.
 205 C
 206       nfl=0
 207       icg=1
 208       do i=1,14
 209         do j=1,14
 210           if (print_order(i).eq.j) then
 211             iw(print_order(i))=j
 212             goto 1121
 213           endif
 214         enddo
 215 1121    continue
 216       enddo
 217       calc_grad=.false.
 218 C Set timers and counters for the respective routines
 219       t_func = 0.0d0
 220       t_grad = 0.0d0
 221       t_fhel = 0.0d0
 222       t_fbet = 0.0d0
 223       t_ghel = 0.0d0
 224       t_gbet = 0.0d0
 225       t_viol = 0.0d0
 226       t_gviol = 0.0d0
 227       n_func = 0
 228       n_grad = 0
 229       n_fhel = 0
 230       n_fbet = 0
 231       n_ghel = 0
 232       n_gbet = 0
 233       n_viol = 0
 234       n_gviol = 0
 235       n_map = 0
 236 #ifndef SPLITELE
 237       nprint_ene=nprint_ene-1
 238 #endif
 239       return
 240       end
 241 c-------------------------------------------------------------------------
 242       block data nazwy
 243       implicit real*8 (a-h,o-z)
 244       include 'DIMENSIONS'
 245       include 'DIMENSIONS.ZSCOPT'
 246       include 'COMMON.NAMES'
 247       include 'COMMON.WEIGHTS'
 248       include 'COMMON.FFIELD'
 249       data restyp /
 250      &'DD' ,'DPR','DLY','DAR','DHI','DAS','DGL','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','D'/
 254       data onelet /
 255      &'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','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=ind_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