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