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