source/cluster/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 'sizesclu.dat'
   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.NAMES"
  22       include "COMMON.TIME1"
  23 C
  24 C The following is just to define auxiliary variables used in angle conversion
  25 C
  26       pi=4.0D0*datan(1.0D0)
  27       dwapi=2.0D0*pi
  28       dwapi3=dwapi/3.0D0
  29       pipol=0.5D0*pi
  30       deg2rad=pi/180.0D0
  31       rad2deg=1.0D0/deg2rad
  32       angmin=10.0D0*deg2rad
  33       Rgas = 1.987D-3
  34 C
  35 C Define I/O units.
  36 C
  37       inp=    1
  38       iout=   2
  39       ipdbin= 3
  40       ipdb=   7
  41       imol2= 18
  42       jplot= 19
  43       jstatin=10
  44       imol2=  4
  45       igeom=  8
  46       intin=  9
  47       ithep= 11
  48       irotam=12
  49       itorp= 13
  50       itordp= 23
  51       ielep= 14
  52       isidep=15
  53       isidep1=22
  54       iscpp=25
  55       icbase=16
  56       ifourier=20
  57       istat= 17
  58       ibond=28
  59       isccor=29
  60       jrms=30
  61 C
  62 C Set default weights of the energy terms.
  63 C
  64       wlong=1.0D0
  65       welec=1.0D0
  66       wtor =1.0D0
  67       wang =1.0D0
  68       wscloc=1.0D0
  69       wstrain=1.0D0
  70 C
  71 C Zero out tables.
  72 C
  73       ndih_constr=0
  74       do i=1,maxres2
  75         do j=1,3
  76           c(j,i)=0.0D0
  77           dc(j,i)=0.0D0
  78         enddo
  79       enddo
  80       do i=1,maxres
  81         do j=1,3
  82           xloc(j,i)=0.0D0
  83         enddo
  84       enddo
  85       do i=1,ntyp
  86         do j=1,ntyp
  87           aa(i,j)=0.0D0
  88           bb(i,j)=0.0D0
  89           augm(i,j)=0.0D0
  90           sigma(i,j)=0.0D0
  91           r0(i,j)=0.0D0
  92           chi(i,j)=0.0D0
  93         enddo
  94         do j=1,2
  95           bad(i,j)=0.0D0
  96         enddo
  97         chip(i)=0.0D0
  98         alp(i)=0.0D0
  99         sigma0(i)=0.0D0
 100         sigii(i)=0.0D0
 101         rr0(i)=0.0D0
 102         a0thet(i)=0.0D0
 103         do j=1,2
 104          do ichir1=-1,1
 105           do ichir2=-1,1
 106           athet(j,i,ichir1,ichir2)=0.0D0
 107           bthet(j,i,ichir1,ichir2)=0.0D0
 108           enddo
 109          enddo
 110         enddo
 111         do j=0,3
 112           polthet(j,i)=0.0D0
 113         enddo
 114         do j=1,3
 115           gthet(j,i)=0.0D0
 116         enddo
 117         theta0(i)=0.0D0
 118         sig0(i)=0.0D0
 119         sigc0(i)=0.0D0
 120         do j=1,maxlob
 121           bsc(j,i)=0.0D0
 122           do k=1,3
 123             censc(k,j,i)=0.0D0
 124           enddo
 125           do k=1,3
 126             do l=1,3
 127               gaussc(l,k,j,i)=0.0D0
 128             enddo
 129           enddo
 130           nlob(i)=0
 131         enddo
 132       enddo
 133       nlob(ntyp1)=0
 134       dsc(ntyp1)=0.0D0
 135       do i=-maxtor,maxtor
 136         itortyp(i)=0
 137        do iblock=1,2
 138         do j=-maxtor,maxtor
 139           do k=1,maxterm
 140             v1(k,j,i,iblock)=0.0D0
 141             v2(k,j,i,iblock)=0.0D0
 142            enddo
 143          enddo
 144        enddo
 145       enddo
 146       do iblock=1,2
 147        do i=-maxtor,maxtor
 148         do j=-maxtor,maxtor
 149          do k=-maxtor,maxtor
 150           do l=1,maxtermd_1
 151             v1c(1,l,i,j,k,iblock)=0.0D0
 152             v1s(1,l,i,j,k,iblock)=0.0D0
 153             v1c(2,l,i,j,k,iblock)=0.0D0
 154             v1s(2,l,i,j,k,iblock)=0.0D0
 155           enddo !l
 156           do l=1,maxtermd_2
 157            do m=1,maxtermd_2
 158             v2c(m,l,i,j,k,iblock)=0.0D0
 159             v2s(m,l,i,j,k,iblock)=0.0D0
 160            enddo !m
 161           enddo !l
 162         enddo !k
 163        enddo !j
 164       enddo !i
 165       enddo !iblock
 166       do i=1,maxres
 167         itype(i)=0
 168         itel(i)=0
 169       enddo
 170 C Initialize the bridge arrays
 171       ns=0
 172       nss=0
 173       nhpb=0
 174       do i=1,maxss
 175         iss(i)=0
 176       enddo
 177       do i=1,maxss
 178         dhpb(i)=0.0D0
 179       enddo
 180       do i=1,maxss
 181         ihpb(i)=0
 182         jhpb(i)=0
 183       enddo
 184 C
 185 C Initialize timing.
 186 C
 187       call set_timers
 188 C
 189 C Initialize variables used in minimization.
 190 C
 191 c     maxfun=5000
 192 c     maxit=2000
 193       maxfun=500
 194       maxit=200
 195       tolf=1.0D-2
 196       rtolf=5.0D-4
 197 C
 198 C Initialize the variables responsible for the mode of gradient storage.
 199 C
 200       nfl=0
 201       icg=1
 202       do i=1,14
 203         do j=1,14
 204           if (print_order(i).eq.j) then
 205             iw(print_order(i))=j
 206             goto 1121
 207           endif
 208         enddo
 209 1121    continue
 210       enddo
 211       calc_grad=.false.
 212 C Set timers and counters for the respective routines
 213       t_func = 0.0d0
 214       t_grad = 0.0d0
 215       t_fhel = 0.0d0
 216       t_fbet = 0.0d0
 217       t_ghel = 0.0d0
 218       t_gbet = 0.0d0
 219       t_viol = 0.0d0
 220       t_gviol = 0.0d0
 221       n_func = 0
 222       n_grad = 0
 223       n_fhel = 0
 224       n_fbet = 0
 225       n_ghel = 0
 226       n_gbet = 0
 227       n_viol = 0
 228       n_gviol = 0
 229       n_map = 0
 230 #ifndef SPLITELE
 231       nprint_ene=nprint_ene-1
 232 #endif
 233       return
 234       end
 235 c-------------------------------------------------------------------------
 236       block data nazwy
 237       implicit real*8 (a-h,o-z)
 238       include 'DIMENSIONS'
 239       include 'sizesclu.dat'
 240       include 'COMMON.NAMES'
 241       include 'COMMON.FFIELD'
 242       data restyp /
 243      &'DD','DAU','DAI','DDB','DSM','DPR','DLY','DAR','DHI','DAS','DGL',
 244      & 'DSG','DGN','DSN','DTH',
 245      &'DYY','DAL','DTY','DTR','DVA','DLE','DIL','DPN','MED','DCY','ZER',
 246      &'CYS','MET','PHE','ILE','LEU','VAL','TRP','TYR','ALA','GLY','THR',
 247      &'SER','GLN','ASN','GLU','ASP','HIS','ARG','LYS','PRO','SME','DBZ',
 248      &'AIB','ABU','D'/
 249       data onelet /
 250      &'z','z','z','z','z','p','k','r','h','d','e','n','q','s','t','g',
 251      &'a','y','w','v','l','i','f','m','c','x',
 252      &'C','M','F','I','L','V','W','Y','A','G','T',
 253      &'S','Q','N','E','D','H','R','K','P','z','z','z','z','X'/
 254       data potname /'LJ','LJK','BP','GB','GBV'/
 255       data ename /
 256      &   "EVDW SC-SC","EVDW2 SC-p","EES p-p","ECORR4 ","ECORR5 ",
 257      &   "ECORR6 ","EELLO ","ETURN3 ","ETURN4 ","ETURN6 ",
 258      &   "EBE bend","ESC SCloc","ETORS ","ETORSD ","EHPB","EVDWPP",
 259      &   "ESTR","ESCCOR","EVDW2_14","","EVDW_T"/
 260       data wname /
 261      &   "WSC","WSCP","WELEC","WCORR","WCORR5","WCORR6","WEL_LOC",
 262      &   "WTURN3","WTURN4","WTURN6","WANG","WSCLOC","WTOR","WTORD",
 263      &   "WHPB","WVDWPP","WBOND","WSCCOR","WSCP14","","WSC"/
 264       data nprint_ene /21/
 265       data print_order /1,2,3,18,11,12,13,14,4,5,6,7,8,9,10,19,
 266      &  16,15,17,20,21/
 267       end
 268 c---------------------------------------------------------------------------
 269       subroutine init_int_table
 270       implicit real*8 (a-h,o-z)
 271       include 'DIMENSIONS'
 272       include 'sizesclu.dat'
 273 #ifdef MPI
 274       include 'mpif.h'
 275 #endif
 276 #ifdef MPL
 277       include 'COMMON.INFO'
 278 #endif
 279       include 'COMMON.CHAIN'
 280       include 'COMMON.INTERACT'
 281       include 'COMMON.LOCAL'
 282       include 'COMMON.SBRIDGE'
 283       include 'COMMON.IOUNITS'
 284       logical scheck,lprint
 285 #ifdef MPL
 286       integer my_sc_int(0:max_fg_Procs-1),my_sc_intt(0:max_fg_Procs),
 287      & my_ele_int(0:max_fg_Procs-1),my_ele_intt(0:max_fg_Procs)
 288 C... Determine the numbers of start and end SC-SC interaction
 289 C... to deal with by current processor.
 290       lprint=.false.
 291       if (lprint)
 292      &write (iout,*) 'INIT_INT_TABLE nres=',nres,' nnt=',nnt,' nct=',nct
 293       n_sc_int_tot=(nct-nnt+1)*(nct-nnt)/2-nss
 294       MyRank=MyID-(MyGroup-1)*fgProcs
 295       call int_bounds(n_sc_int_tot,my_sc_inds,my_sc_inde)
 296       if (lprint)
 297      &  write (iout,*) 'Processor',MyID,' MyRank',MyRank,
 298      &  ' n_sc_int_tot',n_sc_int_tot,' my_sc_inds=',my_sc_inds,
 299      &  ' my_sc_inde',my_sc_inde
 300       ind_sctint=0
 301       iatsc_s=0
 302       iatsc_e=0
 303 #endif
 304       lprint=.false.
 305       do i=1,maxres
 306         nint_gr(i)=0
 307         nscp_gr(i)=0
 308         do j=1,maxint_gr
 309           istart(i,1)=0
 310           iend(i,1)=0
 311           ielstart(i)=0
 312           ielend(i)=0
 313           iscpstart(i,1)=0
 314           iscpend(i,1)=0
 315         enddo
 316       enddo
 317       ind_scint=0
 318       ind_scint_old=0
 319 cd    write (iout,*) 'ns=',ns,' nss=',nss,' ihpb,jhpb',
 320 cd   &   (ihpb(i),jhpb(i),i=1,nss)
 321       do i=nnt,nct-1
 322         scheck=.false.
 323         if (dyn_ss) goto 10
 324         do ii=1,nss
 325           if (ihpb(ii).eq.i+nres) then
 326             scheck=.true.
 327             jj=jhpb(ii)-nres
 328             goto 10
 329           endif
 330         enddo
 331    10   continue
 332 cd      write (iout,*) 'i=',i,' scheck=',scheck,' jj=',jj
 333         if (scheck) then
 334           if (jj.eq.i+1) then
 335 #ifdef MPL
 336             write (iout,*) 'jj=i+1'
 337             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
 338      & iatsc_s,iatsc_e,i+2,nct,nint_gr(i),istart(i,1),iend(i,1),*12)
 339 #else
 340             nint_gr(i)=1
 341             istart(i,1)=i+2
 342             iend(i,1)=nct
 343 #endif
 344           else if (jj.eq.nct) then
 345 #ifdef MPL
 346             write (iout,*) 'jj=nct'
 347             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
 348      &  iatsc_s,iatsc_e,i+1,nct-1,nint_gr(i),istart(i,1),iend(i,1),*12)
 349 #else
 350             nint_gr(i)=1
 351             istart(i,1)=i+1
 352             iend(i,1)=nct-1
 353 #endif
 354           else
 355 #ifdef MPL
 356             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
 357      & iatsc_s,iatsc_e,i+1,jj-1,nint_gr(i),istart(i,1),iend(i,1),*12)
 358             ii=nint_gr(i)+1
 359             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
 360      & iatsc_s,iatsc_e,jj+1,nct,nint_gr(i),istart(i,ii),iend(i,ii),*12)
 361 #else
 362             nint_gr(i)=2
 363             istart(i,1)=i+1
 364             iend(i,1)=jj-1
 365             istart(i,2)=jj+1
 366             iend(i,2)=nct
 367 #endif
 368           endif
 369         else
 370 #ifdef MPL
 371           call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
 372      &    iatsc_s,iatsc_e,i+1,nct,nint_gr(i),istart(i,1),iend(i,1),*12)
 373 #else
 374           nint_gr(i)=1
 375           istart(i,1)=i+1
 376           iend(i,1)=nct
 377           ind_scint=int_scint+nct-i
 378 #endif
 379         endif
 380 #ifdef MPL
 381         ind_scint_old=ind_scint
 382 #endif
 383       enddo
 384    12 continue
 385 #ifndef MPL
 386       iatsc_s=nnt
 387       iatsc_e=nct-1
 388 #endif
 389 #ifdef MPL
 390       if (lprint) then
 391         write (iout,*) 'Processor',MyID,' Group',MyGroup
 392         write (iout,*) 'iatsc_s=',iatsc_s,' iatsc_e=',iatsc_e
 393       endif
 394 #endif
 395       if (lprint) then
 396       write (iout,'(a)') 'Interaction array:'
 397       do i=iatsc_s,iatsc_e
 398         write (iout,'(i3,2(2x,2i3))')
 399      & i,(istart(i,iint),iend(i,iint),iint=1,nint_gr(i))
 400       enddo
 401       endif
 402       ispp=2
 403 #ifdef MPL
 404 C Now partition the electrostatic-interaction array
 405       npept=nct-nnt
 406       nele_int_tot=(npept-ispp)*(npept-ispp+1)/2
 407       call int_bounds(nele_int_tot,my_ele_inds,my_ele_inde)
 408       if (lprint)
 409      & write (iout,*) 'Processor',MyID,' MyRank',MyRank,
 410      &  ' nele_int_tot',nele_int_tot,' my_ele_inds=',my_ele_inds,
 411      &               ' my_ele_inde',my_ele_inde
 412       iatel_s=0
 413       iatel_e=0
 414       ind_eleint=0
 415       ind_eleint_old=0
 416       do i=nnt,nct-3
 417         ijunk=0
 418         call int_partition(ind_eleint,my_ele_inds,my_ele_inde,i,
 419      &    iatel_s,iatel_e,i+ispp,nct-1,ijunk,ielstart(i),ielend(i),*13)
 420       enddo ! i
 421    13 continue
 422 #else
 423       iatel_s=nnt
 424       iatel_e=nct-3
 425       do i=iatel_s,iatel_e
 426         ielstart(i)=i+2
 427         ielend(i)=nct-1
 428       enddo
 429 #endif
 430       if (lprint) then
 431         write (iout,'(a)') 'Electrostatic interaction array:'
 432         do i=iatel_s,iatel_e
 433           write (iout,'(i3,2(2x,2i3))') i,ielstart(i),ielend(i)
 434         enddo
 435       endif ! lprint
 436 c     iscp=3
 437       iscp=2
 438 C Partition the SC-p interaction array
 439 #ifdef MPL
 440       nscp_int_tot=(npept-iscp+1)*(npept-iscp+1)
 441       call int_bounds(nscp_int_tot,my_scp_inds,my_scp_inde)
 442       if (lprint)
 443      & write (iout,*) 'Processor',MyID,' MyRank',MyRank,
 444      &  ' nscp_int_tot',nscp_int_tot,' my_scp_inds=',my_scp_inds,
 445      &               ' my_scp_inde',my_scp_inde
 446       iatscp_s=0
 447       iatscp_e=0
 448       ind_scpint=0
 449       ind_scpint_old=0
 450       do i=nnt,nct-1
 451         if (i.lt.nnt+iscp) then
 452 cd        write (iout,*) 'i.le.nnt+iscp'
 453           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
 454      &      iatscp_s,iatscp_e,i+iscp,nct,nscp_gr(i),iscpstart(i,1),
 455      &      iscpend(i,1),*14)
 456         else if (i.gt.nct-iscp) then
 457 cd        write (iout,*) 'i.gt.nct-iscp'
 458           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
 459      &      iatscp_s,iatscp_e,nnt,i-iscp,nscp_gr(i),iscpstart(i,1),
 460      &      iscpend(i,1),*14)
 461         else
 462           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
 463      &      iatscp_s,iatscp_e,nnt,i-iscp,nscp_gr(i),iscpstart(i,1),
 464      &      iscpend(i,1),*14)
 465           ii=nscp_gr(i)+1
 466           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
 467      &      iatscp_s,iatscp_e,i+iscp,nct,nscp_gr(i),iscpstart(i,ii),
 468      &      iscpend(i,ii),*14)
 469         endif
 470       enddo ! i
 471    14 continue
 472 #else
 473       iatscp_s=nnt
 474       iatscp_e=nct-1
 475       do i=nnt,nct-1
 476         if (i.lt.nnt+iscp) then
 477           nscp_gr(i)=1
 478           iscpstart(i,1)=i+iscp
 479           iscpend(i,1)=nct
 480         elseif (i.gt.nct-iscp) then
 481           nscp_gr(i)=1
 482           iscpstart(i,1)=nnt
 483           iscpend(i,1)=i-iscp
 484         else
 485           nscp_gr(i)=2
 486           iscpstart(i,1)=nnt
 487           iscpend(i,1)=i-iscp
 488           iscpstart(i,2)=i+iscp
 489           iscpend(i,2)=nct
 490         endif
 491       enddo ! i
 492 #endif
 493       if (lprint) then
 494         write (iout,'(a)') 'SC-p interaction array:'
 495         do i=iatscp_s,iatscp_e
 496           write (iout,'(i3,2(2x,2i3))')
 497      &         i,(iscpstart(i,j),iscpend(i,j),j=1,nscp_gr(i))
 498         enddo
 499       endif ! lprint
 500 C Partition local interactions
 501 #ifdef MPL
 502       call int_bounds(nres-2,loc_start,loc_end)
 503       loc_start=loc_start+1
 504       loc_end=loc_end+1
 505       call int_bounds(nres-2,ithet_start,ithet_end)
 506       ithet_start=ithet_start+2
 507       ithet_end=ithet_end+2
 508       call int_bounds(nct-nnt-2,iphi_start,iphi_end)
 509       iphi_start=iphi_start+nnt+2
 510       iphi_end=iphi_end+nnt+2
 511       call int_bounds(nres-3,itau_start,itau_end)
 512       itau_start=itau_start+3
 513       itau_end=itau_end+3
 514       if (lprint) then
 515         write (iout,*) 'Processor:',MyID,
 516      & ' loc_start',loc_start,' loc_end',loc_end,
 517      & ' ithet_start',ithet_start,' ithet_end',ithet_end,
 518      & ' iphi_start',iphi_start,' iphi_end',iphi_end
 519         write (*,*) 'Processor:',MyID,
 520      & ' loc_start',loc_start,' loc_end',loc_end,
 521      & ' ithet_start',ithet_start,' ithet_end',ithet_end,
 522      & ' iphi_start',iphi_start,' iphi_end',iphi_end
 523       endif
 524       if (fgprocs.gt.1 .and. MyID.eq.BossID) then
 525         write(iout,'(i10,a,i10,a,i10,a/a,i3,a)') n_sc_int_tot,' SC-SC ',
 526      & nele_int_tot,' electrostatic and ',nscp_int_tot,
 527      & ' SC-p interactions','were distributed among',fgprocs,
 528      & ' fine-grain processors.'
 529       endif
 530 #else
 531       loc_start=2
 532       loc_end=nres-1
 533       ithet_start=3
 534       ithet_end=nres
 535       iphi_start=nnt+3
 536       iphi_end=nct
 537       itau_start=4
 538       itau_end=nres
 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