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