source/unres/src_MD-M/readrtns_CSA.F

   1       subroutine readrtns
   2       implicit real*8 (a-h,o-z)
   3       include 'DIMENSIONS'
   4 #ifdef MPI
   5       include 'mpif.h'
   6 #endif
   7       include 'COMMON.SETUP'
   8       include 'COMMON.CONTROL'
   9       include 'COMMON.SBRIDGE'
  10       include 'COMMON.IOUNITS'
  11       include 'COMMON.SPLITELE'
  12       logical file_exist
  13 C Read force-field parameters except weights
  14       call parmread
  15 C Read job setup parameters
  16       call read_control
  17 C Read control parameters for energy minimzation if required
  18       if (minim) call read_minim
  19 C Read MCM control parameters if required
  20       if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
  21 C Read MD control parameters if reqjuired
  22       if (modecalc.eq.12) call read_MDpar
  23 C Read MREMD control parameters if required
  24       if (modecalc.eq.14) then
  25          call read_MDpar
  26          call read_REMDpar
  27       endif
  28 C Read MUCA control parameters if required
  29       if (lmuca) call read_muca
  30 C Read CSA control parameters if required (from fort.40 if exists
  31 C otherwise from general input file)
  32       if (modecalc.eq.8) then
  33        inquire (file="fort.40",exist=file_exist)
  34        if (.not.file_exist) call csaread
  35       endif
  36 cfmc      if (modecalc.eq.10) call mcmfread
  37 C Read molecule information, molecule geometry, energy-term weights, and
  38 C restraints if requested
  39       call molread
  40 C Print restraint information
  41 #ifdef MPI
  42       if (.not. out1file .or. me.eq.king) then
  43 #endif
  44       if (nhpb.gt.nss)
  45      &write (iout,'(a,i5,a)') "The following",nhpb-nss,
  46      & " distance constraints have been imposed"
  47       do i=nss+1,nhpb
  48         write (iout,'(3i6,f10.5)') i-nss,ihpb(i),jhpb(i),forcon(i)
  49       enddo
  50 #ifdef MPI
  51       endif
  52 #endif
  53 c      print *,"Processor",myrank," leaves READRTNS"
  54       return
  55       end
  56 C-------------------------------------------------------------------------------
  57       subroutine read_control
  58 C
  59 C Read contorl data
  60 C
  61       implicit real*8 (a-h,o-z)
  62       include 'DIMENSIONS'
  63 #ifdef MP
  64       include 'mpif.h'
  65       logical OKRandom, prng_restart
  66       real*8  r1
  67 #endif
  68       include 'COMMON.IOUNITS'
  69       include 'COMMON.TIME1'
  70       include 'COMMON.THREAD'
  71       include 'COMMON.SBRIDGE'
  72       include 'COMMON.CONTROL'
  73       include 'COMMON.MCM'
  74       include 'COMMON.MAP'
  75       include 'COMMON.HEADER'
  76       include 'COMMON.CSA'
  77       include 'COMMON.CHAIN'
  78       include 'COMMON.MUCA'
  79       include 'COMMON.MD'
  80       include 'COMMON.FFIELD'
  81       include 'COMMON.INTERACT'
  82       include 'COMMON.SETUP'
  83       include 'COMMON.SPLITELE'
  84       include 'COMMON.SHIELD'
  85       COMMON /MACHSW/ KDIAG,ICORFL,IXDR
  86       character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
  87       character*80 ucase
  88       character*320 controlcard
  89
  90       nglob_csa=0
  91       eglob_csa=1d99
  92       nmin_csa=0
  93       read (INP,'(a)') titel
  94       call card_concat(controlcard)
  95 c      out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
  96 c      print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
  97       call reada(controlcard,'SEED',seed,0.0D0)
  98       call random_init(seed)
  99 C Set up the time limit (caution! The time must be input in minutes!)
 100       read_cart=index(controlcard,'READ_CART').gt.0
 101       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
 102 C this variable with_theta_constr is the variable which allow to read and execute the
 103 C constrains on theta angles WITH_THETA_CONSTR is the keyword
 104       with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
 105       write (iout,*) "with_theta_constr ",with_theta_constr
 106       call readi(controlcard,'SYM',symetr,1)
 107       call reada(controlcard,'TIMLIM',timlim,2800.0D0) ! default 16 hours
 108       unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
 109       call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
 110       call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
 111       call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
 112       call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
 113       call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
 114       call reada(controlcard,'DRMS',drms,0.1D0)
 115       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 116        write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
 117        write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
 118        write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
 119        write (iout,'(a,f10.1)')'DRMS    = ',drms
 120        write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
 121        write (iout,'(a,f10.1)') 'Time limit (min):',timlim
 122       endif
 123       call readi(controlcard,'NZ_START',nz_start,0)
 124       call readi(controlcard,'NZ_END',nz_end,0)
 125       call readi(controlcard,'IZ_SC',iz_sc,0)
 126       timlim=60.0D0*timlim
 127       safety = 60.0d0*safety
 128       timem=timlim
 129       modecalc=0
 130       call reada(controlcard,"T_BATH",t_bath,300.0d0)
 131       minim=(index(controlcard,'MINIMIZE').gt.0)
 132       dccart=(index(controlcard,'CART').gt.0)
 133       overlapsc=(index(controlcard,'OVERLAP').gt.0)
 134       overlapsc=.not.overlapsc
 135       searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
 136       searchsc=.not.searchsc
 137       sideadd=(index(controlcard,'SIDEADD').gt.0)
 138       energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
 139       outpdb=(index(controlcard,'PDBOUT').gt.0)
 140       outmol2=(index(controlcard,'MOL2OUT').gt.0)
 141       pdbref=(index(controlcard,'PDBREF').gt.0)
 142       refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
 143       indpdb=index(controlcard,'PDBSTART')
 144       extconf=(index(controlcard,'EXTCONF').gt.0)
 145       AFMlog=(index(controlcard,'AFM'))
 146       selfguide=(index(controlcard,'SELFGUIDE'))
 147       print *,'AFMlog',AFMlog,selfguide,"KUPA"
 148       call readi(controlcard,'TUBEMOD',tubelog,0)
 149       write (iout,*) TUBElog,"TUBEMODE"
 150       call readi(controlcard,'GENCONSTR',genconstr,0)
 151 C      write (iout,*) TUBElog,"TUBEMODE"
 152       call readi(controlcard,'IPRINT',iprint,0)
 153 C SHIELD keyword sets if the shielding effect of side-chains is used
 154 C 0 denots no shielding is used all peptide are equally despite the
 155 C solvent accesible area
 156 C 1 the newly introduced function
 157 C 2 reseved for further possible developement
 158       call readi(controlcard,'SHIELD',shield_mode,0)
 159 C      if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 160         write(iout,*) "shield_mode",shield_mode
 161 C      endif
 162       call readi(controlcard,'TORMODE',tor_mode,0)
 163 C      if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 164         write(iout,*) "torsional and valence angle mode",tor_mode
 165       call readi(controlcard,'MAXGEN',maxgen,10000)
 166       call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
 167       call readi(controlcard,"KDIAG",kdiag,0)
 168       call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
 169       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
 170      & write (iout,*) "RESCALE_MODE",rescale_mode
 171       split_ene=index(controlcard,'SPLIT_ENE').gt.0
 172       if (index(controlcard,'REGULAR').gt.0.0D0) then
 173         call reada(controlcard,'WEIDIS',weidis,0.1D0)
 174         modecalc=1
 175         refstr=.true.
 176       endif
 177       if (index(controlcard,'CHECKGRAD').gt.0) then
 178         modecalc=5
 179         if (index(controlcard,'CART').gt.0) then
 180           icheckgrad=1
 181         elseif (index(controlcard,'CARINT').gt.0) then
 182           icheckgrad=2
 183         else
 184           icheckgrad=3
 185         endif
 186         call reada(controlcard,'DELTA',aincr,1.0d-7)
 187       elseif (index(controlcard,'THREAD').gt.0) then
 188         modecalc=2
 189         call readi(controlcard,'THREAD',nthread,0)
 190         if (nthread.gt.0) then
 191           call reada(controlcard,'WEIDIS',weidis,0.1D0)
 192         else
 193           if (fg_rank.eq.0)
 194      &    write (iout,'(a)')'A number has to follow the THREAD keyword.'
 195           stop 'Error termination in Read_Control.'
 196         endif
 197       else if (index(controlcard,'MCMA').gt.0) then
 198         modecalc=3
 199       else if (index(controlcard,'MCEE').gt.0) then
 200         modecalc=6
 201       else if (index(controlcard,'MULTCONF').gt.0) then
 202         modecalc=4
 203       else if (index(controlcard,'MAP').gt.0) then
 204         modecalc=7
 205         call readi(controlcard,'MAP',nmap,0)
 206       else if (index(controlcard,'CSA').gt.0) then
 207         modecalc=8
 208 crc      else if (index(controlcard,'ZSCORE').gt.0) then
 209 crc
 210 crc  ZSCORE is rm from UNRES, modecalc=9 is available
 211 crc
 212 crc        modecalc=9
 213 cfcm      else if (index(controlcard,'MCMF').gt.0) then
 214 cfmc        modecalc=10
 215       else if (index(controlcard,'SOFTREG').gt.0) then
 216         modecalc=11
 217       else if (index(controlcard,'CHECK_BOND').gt.0) then
 218         modecalc=-1
 219       else if (index(controlcard,'TEST').gt.0) then
 220         modecalc=-2
 221       else if (index(controlcard,'MD').gt.0) then
 222         modecalc=12
 223       else if (index(controlcard,'RE ').gt.0) then
 224         modecalc=14
 225       endif
 226
 227       lmuca=index(controlcard,'MUCA').gt.0
 228       call readi(controlcard,'MUCADYN',mucadyn,0)
 229       call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
 230       if (lmuca .and. (me.eq.king .or. .not.out1file ))
 231      & then
 232        write (iout,*) 'MUCADYN=',mucadyn
 233        write (iout,*) 'MUCASMOOTH=',muca_smooth
 234       endif
 235
 236       iscode=index(controlcard,'ONE_LETTER')
 237       indphi=index(controlcard,'PHI')
 238       indback=index(controlcard,'BACK')
 239       iranconf=index(controlcard,'RAND_CONF')
 240       i2ndstr=index(controlcard,'USE_SEC_PRED')
 241       gradout=index(controlcard,'GRADOUT').gt.0
 242       gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
 243 C  DISTCHAINMAX become obsolete for periodic boundry condition
 244       call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
 245 C Reading the dimensions of box in x,y,z coordinates
 246       call reada(controlcard,'BOXX',boxxsize,100.0d0)
 247       call reada(controlcard,'BOXY',boxysize,100.0d0)
 248       call reada(controlcard,'BOXZ',boxzsize,100.0d0)
 249 c Cutoff range for interactions
 250       call reada(controlcard,"R_CUT",r_cut,15.0d0)
 251       call reada(controlcard,"LAMBDA",rlamb,0.3d0)
 252       call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
 253       call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
 254       if (lipthick.gt.0.0d0) then
 255        bordliptop=(boxzsize+lipthick)/2.0
 256        bordlipbot=bordliptop-lipthick
 257 C      endif
 258       if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
 259      & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
 260       buflipbot=bordlipbot+lipbufthick
 261       bufliptop=bordliptop-lipbufthick
 262       if ((lipbufthick*2.0d0).gt.lipthick)
 263      &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
 264       endif
 265       write(iout,*) "bordliptop=",bordliptop
 266       write(iout,*) "bordlipbot=",bordlipbot
 267       write(iout,*) "bufliptop=",bufliptop
 268       write(iout,*) "buflipbot=",buflipbot
 269       write (iout,*) "SHIELD MODE",shield_mode
 270       if (TUBElog.gt.0) then
 271        call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
 272        call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
 273        call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0)
 274        call reada(controlcard,"RTUBE",tubeR0,0.0d0)
 275        call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
 276        call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
 277        call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
 278        buftubebot=bordtubebot+tubebufthick
 279        buftubetop=bordtubetop-tubebufthick
 280       endif
 281       if (shield_mode.gt.0) then
 282       pi=3.141592d0
 283 C VSolvSphere the volume of solving sphere
 284 C      print *,pi,"pi"
 285 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
 286 C there will be no distinction between proline peptide group and normal peptide
 287 C group in case of shielding parameters
 288       VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
 289       VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
 290       write (iout,*) VSolvSphere,VSolvSphere_div
 291 C long axis of side chain
 292       do i=1,ntyp
 293       long_r_sidechain(i)=vbldsc0(1,i)
 294       short_r_sidechain(i)=sigma0(i)
 295       enddo
 296       buff_shield=1.0d0
 297       endif
 298       if (me.eq.king .or. .not.out1file )
 299      &  write (iout,*) "DISTCHAINMAX",distchainmax
 300
 301       if(me.eq.king.or..not.out1file)
 302      & write (iout,'(2a)') diagmeth(kdiag),
 303      &  ' routine used to diagonalize matrices.'
 304       return
 305       end
 306 c--------------------------------------------------------------------------
 307       subroutine read_REMDpar
 308 C
 309 C Read REMD settings
 310 C
 311       implicit real*8 (a-h,o-z)
 312       include 'DIMENSIONS'
 313       include 'COMMON.IOUNITS'
 314       include 'COMMON.TIME1'
 315       include 'COMMON.MD'
 316 #ifndef LANG0
 317       include 'COMMON.LANGEVIN'
 318 #else
 319       include 'COMMON.LANGEVIN.lang0'
 320 #endif
 321       include 'COMMON.INTERACT'
 322       include 'COMMON.NAMES'
 323       include 'COMMON.GEO'
 324       include 'COMMON.REMD'
 325       include 'COMMON.CONTROL'
 326       include 'COMMON.SETUP'
 327       character*80 ucase
 328       character*320 controlcard
 329       character*3200 controlcard1
 330       integer iremd_m_total
 331
 332       if(me.eq.king.or..not.out1file)
 333      & write (iout,*) "REMD setup"
 334
 335       call card_concat(controlcard)
 336       call readi(controlcard,"NREP",nrep,3)
 337       call readi(controlcard,"NSTEX",nstex,1000)
 338       call reada(controlcard,"RETMIN",retmin,10.0d0)
 339       call reada(controlcard,"RETMAX",retmax,1000.0d0)
 340       mremdsync=(index(controlcard,'SYNC').gt.0)
 341       call readi(controlcard,"NSYN",i_sync_step,100)
 342       restart1file=(index(controlcard,'REST1FILE').gt.0)
 343       traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
 344       call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
 345       if(max_cache_traj_use.gt.max_cache_traj)
 346      &           max_cache_traj_use=max_cache_traj
 347       if(me.eq.king.or..not.out1file) then
 348 cd       if (traj1file) then
 349 crc caching is in testing - NTWX is not ignored
 350 cd        write (iout,*) "NTWX value is ignored"
 351 cd        write (iout,*) "  trajectory is stored to one file by master"
 352 cd        write (iout,*) "  before exchange at NSTEX intervals"
 353 cd       endif
 354        write (iout,*) "NREP= ",nrep
 355        write (iout,*) "NSTEX= ",nstex
 356        write (iout,*) "SYNC= ",mremdsync
 357        write (iout,*) "NSYN= ",i_sync_step
 358        write (iout,*) "TRAJCACHE= ",max_cache_traj_use
 359       endif
 360       remd_tlist=.false.
 361       if (index(controlcard,'TLIST').gt.0) then
 362          remd_tlist=.true.
 363          call card_concat(controlcard1)
 364          read(controlcard1,*) (remd_t(i),i=1,nrep)
 365          if(me.eq.king.or..not.out1file)
 366      &    write (iout,*)'tlist',(remd_t(i),i=1,nrep)
 367       endif
 368       remd_mlist=.false.
 369       if (index(controlcard,'MLIST').gt.0) then
 370          remd_mlist=.true.
 371          call card_concat(controlcard1)
 372          read(controlcard1,*) (remd_m(i),i=1,nrep)
 373          if(me.eq.king.or..not.out1file) then
 374           write (iout,*)'mlist',(remd_m(i),i=1,nrep)
 375           iremd_m_total=0
 376           do i=1,nrep
 377            iremd_m_total=iremd_m_total+remd_m(i)
 378           enddo
 379            write (iout,*) 'Total number of replicas ',iremd_m_total
 380           endif
 381          endif
 382       if(me.eq.king.or..not.out1file)
 383      &   write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
 384       return
 385       end
 386 c--------------------------------------------------------------------------
 387       subroutine read_MDpar
 388 C
 389 C Read MD settings
 390 C
 391       implicit real*8 (a-h,o-z)
 392       include 'DIMENSIONS'
 393       include 'COMMON.IOUNITS'
 394       include 'COMMON.TIME1'
 395       include 'COMMON.MD'
 396 #ifndef LANG0
 397       include 'COMMON.LANGEVIN'
 398 #else
 399       include 'COMMON.LANGEVIN.lang0'
 400 #endif
 401       include 'COMMON.INTERACT'
 402       include 'COMMON.NAMES'
 403       include 'COMMON.GEO'
 404       include 'COMMON.SETUP'
 405       include 'COMMON.CONTROL'
 406       include 'COMMON.SPLITELE'
 407       character*80 ucase
 408       character*320 controlcard
 409
 410       call card_concat(controlcard)
 411       call readi(controlcard,"NSTEP",n_timestep,1000000)
 412       call readi(controlcard,"NTWE",ntwe,100)
 413       call readi(controlcard,"NTWX",ntwx,1000)
 414       call reada(controlcard,"DT",d_time,1.0d-1)
 415       call reada(controlcard,"DVMAX",dvmax,2.0d1)
 416       call reada(controlcard,"DAMAX",damax,1.0d1)
 417       call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
 418       call readi(controlcard,"LANG",lang,0)
 419       RESPA = index(controlcard,"RESPA") .gt. 0
 420       call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
 421       ntime_split0=ntime_split
 422       call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
 423       ntime_split0=ntime_split
 424 c      call reada(controlcard,"R_CUT",r_cut,2.0d0)
 425 c      call reada(controlcard,"LAMBDA",rlamb,0.3d0)
 426       rest = index(controlcard,"REST").gt.0
 427       tbf = index(controlcard,"TBF").gt.0
 428       tnp = index(controlcard,"NOSEPOINCARE99").gt.0
 429       tnp1 = index(controlcard,"NOSEPOINCARE01").gt.0
 430       tnh = index(controlcard,"NOSEHOOVER96").gt.0
 431       if (RESPA.and.tnh)then
 432         xiresp = index(controlcard,"XIRESP").gt.0
 433       endif
 434       write(iout,*) "tnh",tnh
 435       call reada(controlcard,"Q_NP",Q_np,0.1d0)
 436       usampl = index(controlcard,"USAMPL").gt.0
 437
 438       mdpdb = index(controlcard,"MDPDB").gt.0
 439       call reada(controlcard,"T_BATH",t_bath,300.0d0)
 440       call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
 441       call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
 442       call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
 443       if (count_reset_moment.eq.0) count_reset_moment=1000000000
 444       call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
 445       reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
 446       reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
 447       if (count_reset_vel.eq.0) count_reset_vel=1000000000
 448       large = index(controlcard,"LARGE").gt.0
 449       print_compon = index(controlcard,"PRINT_COMPON").gt.0
 450       rattle = index(controlcard,"RATTLE").gt.0
 451 c  if performing umbrella sampling, fragments constrained are read from the fragment file
 452       nset=0
 453       if(usampl) then
 454         call read_fragments
 455       endif
 456
 457       if(me.eq.king.or..not.out1file) then
 458        write (iout,*)
 459        write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
 460        write (iout,*)
 461        write (iout,'(a)') "The units are:"
 462        write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
 463        write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
 464      &  " acceleration: angstrom/(48.9 fs)**2"
 465        write (iout,'(a)') "energy: kcal/mol, temperature: K"
 466        write (iout,*)
 467        write (iout,'(a60,i10)') "Number of time steps:",n_timestep
 468        write (iout,'(a60,f10.5,a)')
 469      &  "Initial time step of numerical integration:",d_time,
 470      &  " natural units"
 471        write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
 472        if (RESPA) then
 473         write (iout,'(2a,i4,a)')
 474      &    "A-MTS algorithm used; initial time step for fast-varying",
 475      &    " short-range forces split into",ntime_split," steps."
 476         write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
 477      &   r_cut," lambda",rlamb
 478        endif
 479        write (iout,'(2a,f10.5)')
 480      &  "Maximum acceleration threshold to reduce the time step",
 481      &  "/increase split number:",damax
 482        write (iout,'(2a,f10.5)')
 483      &  "Maximum predicted energy drift to reduce the timestep",
 484      &  "/increase split number:",edriftmax
 485        write (iout,'(a60,f10.5)')
 486      & "Maximum velocity threshold to reduce velocities:",dvmax
 487        write (iout,'(a60,i10)') "Frequency of property output:",ntwe
 488        write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
 489        if (rattle) write (iout,'(a60)')
 490      &  "Rattle algorithm used to constrain the virtual bonds"
 491       endif
 492       reset_fricmat=1000
 493       if (lang.gt.0) then
 494         call reada(controlcard,"ETAWAT",etawat,0.8904d0)
 495         call reada(controlcard,"RWAT",rwat,1.4d0)
 496         call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
 497         surfarea=index(controlcard,"SURFAREA").gt.0
 498         call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
 499         if(me.eq.king.or..not.out1file)then
 500          write (iout,'(/a,$)') "Langevin dynamics calculation"
 501          if (lang.eq.1) then
 502           write (iout,'(a/)')
 503      &      " with direct integration of Langevin equations"
 504          else if (lang.eq.2) then
 505           write (iout,'(a/)') " with TINKER stochasic MD integrator"
 506          else if (lang.eq.3) then
 507           write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
 508          else if (lang.eq.4) then
 509           write (iout,'(a/)') " in overdamped mode"
 510          else
 511           write (iout,'(//a,i5)')
 512      &      "=========== ERROR: Unknown Langevin dynamics mode:",lang
 513           stop
 514          endif
 515          write (iout,'(a60,f10.5)') "Temperature:",t_bath
 516          write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
 517          write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
 518          write (iout,'(a60,f10.5)')
 519      &   "Scaling factor of the friction forces:",scal_fric
 520          if (surfarea) write (iout,'(2a,i10,a)')
 521      &     "Friction coefficients will be scaled by solvent-accessible",
 522      &     " surface area every",reset_fricmat," steps."
 523         endif
 524 c Calculate friction coefficients and bounds of stochastic forces
 525         eta=6*pi*cPoise*etawat
 526         if(me.eq.king.or..not.out1file)
 527      &   write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
 528      &   ,eta
 529         gamp=scal_fric*(pstok+rwat)*eta
 530         stdfp=dsqrt(2*Rb*t_bath/d_time)
 531         do i=1,ntyp
 532           gamsc(i)=scal_fric*(restok(i)+rwat)*eta
 533           stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
 534         enddo
 535         if(me.eq.king.or..not.out1file)then
 536          write (iout,'(/2a/)')
 537      &   "Radii of site types and friction coefficients and std's of",
 538      &   " stochastic forces of fully exposed sites"
 539          write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
 540          do i=1,ntyp
 541           write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
 542      &     gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
 543          enddo
 544         endif
 545       else if (tbf) then
 546         if(me.eq.king.or..not.out1file)then
 547          write (iout,'(a)') "Berendsen bath calculation"
 548          write (iout,'(a60,f10.5)') "Temperature:",t_bath
 549          write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
 550          if (reset_moment)
 551      &   write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
 552      &   count_reset_moment," steps"
 553          if (reset_vel)
 554      &    write (iout,'(a,i10,a)')
 555      &    "Velocities will be reset at random every",count_reset_vel,
 556      &   " steps"
 557         endif
 558       else if (tnp .or. tnp1 .or. tnh) then
 559         if (tnp .or. tnp1) then
 560            write (iout,'(a)') "Nose-Poincare bath calculation"
 561            if (tnp) write (iout,'(a)')
 562      & "J.Comput.Phys. 151 114 (1999) S.D.Bond B.J.Leimkuhler B.B.Laird"
 563            if (tnp1) write (iout,'(a)') "JPSJ 70 75 (2001) S. Nose"
 564         else
 565            write (iout,'(a)') "Nose-Hoover bath calculation"
 566            write (iout,'(a)') "Mol.Phys. 87 1117 (1996) Martyna et al."
 567               nresn=1
 568               nyosh=1
 569               nnos=1
 570               do i=1,nnos
 571                qmass(i)=Q_np
 572                xlogs(i)=1.0
 573                vlogs(i)=0.0
 574               enddo
 575               do i=1,nyosh
 576                WDTI(i) = 1.0*d_time/nresn
 577                WDTI2(i)=WDTI(i)/2
 578                WDTI4(i)=WDTI(i)/4
 579                WDTI8(i)=WDTI(i)/8
 580               enddo
 581               if (RESPA) then
 582                if(xiresp) then
 583                  write (iout,'(a)') "NVT-XI-RESPA algorithm"
 584                else
 585                  write (iout,'(a)') "NVT-XO-RESPA algorithm"
 586                endif
 587                do i=1,nyosh
 588                 WDTIi(i) = 1.0*d_time/nresn/ntime_split
 589                 WDTIi2(i)=WDTIi(i)/2
 590                 WDTIi4(i)=WDTIi(i)/4
 591                 WDTIi8(i)=WDTIi(i)/8
 592                enddo
 593               endif
 594         endif
 595       else
 596         if(me.eq.king.or..not.out1file)
 597      &   write (iout,'(a31)') "Microcanonical mode calculation"
 598       endif
 599       if(me.eq.king.or..not.out1file)then
 600        if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
 601        if (usampl) then
 602           write(iout,*) "MD running with constraints."
 603           write(iout,*) "Equilibration time ", eq_time, " mtus."
 604           write(iout,*) "Constraining ", nfrag," fragments."
 605           write(iout,*) "Length of each fragment, weight and q0:"
 606           do iset=1,nset
 607            write (iout,*) "Set of restraints #",iset
 608            do i=1,nfrag
 609               write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
 610      &           ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
 611            enddo
 612            write(iout,*) "constraints between ", npair, "fragments."
 613            write(iout,*) "constraint pairs, weights and q0:"
 614            do i=1,npair
 615             write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
 616      &             ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
 617            enddo
 618            write(iout,*) "angle constraints within ", nfrag_back,
 619      &      "backbone fragments."
 620            write(iout,*) "fragment, weights:"
 621            do i=1,nfrag_back
 622             write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
 623      &         ifrag_back(2,i,iset),wfrag_back(1,i,iset),
 624      &         wfrag_back(2,i,iset),wfrag_back(3,i,iset)
 625            enddo
 626           enddo
 627         iset=mod(kolor,nset)+1
 628        endif
 629       endif
 630       if(me.eq.king.or..not.out1file)
 631      & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
 632       return
 633       end
 634 c------------------------------------------------------------------------------
 635       subroutine molread
 636 C
 637 C Read molecular data.
 638 C
 639       implicit real*8 (a-h,o-z)
 640       include 'DIMENSIONS'
 641 #ifdef MPI
 642       include 'mpif.h'
 643       integer error_msg
 644 #endif
 645       include 'COMMON.IOUNITS'
 646       include 'COMMON.GEO'
 647       include 'COMMON.VAR'
 648       include 'COMMON.INTERACT'
 649       include 'COMMON.LOCAL'
 650       include 'COMMON.NAMES'
 651       include 'COMMON.CHAIN'
 652       include 'COMMON.FFIELD'
 653       include 'COMMON.SBRIDGE'
 654       include 'COMMON.HEADER'
 655       include 'COMMON.CONTROL'
 656       include 'COMMON.DBASE'
 657       include 'COMMON.THREAD'
 658       include 'COMMON.CONTACTS'
 659       include 'COMMON.TORCNSTR'
 660       include 'COMMON.TIME1'
 661       include 'COMMON.BOUNDS'
 662       include 'COMMON.MD'
 663       include 'COMMON.SETUP'
 664       include 'COMMON.SHIELD'
 665       character*4 sequence(maxres)
 666       integer rescode
 667       double precision x(maxvar)
 668       character*256 pdbfile
 669       character*400 weightcard
 670       character*80 weightcard_t,ucase
 671       dimension itype_pdb(maxres)
 672       common /pizda/ itype_pdb
 673       logical seq_comp,fail
 674       double precision energia(0:n_ene)
 675       integer ilen
 676       external ilen
 677 C
 678 C Body
 679 C
 680 C Read weights of the subsequent energy terms.
 681        call card_concat(weightcard)
 682        call reada(weightcard,'WLONG',wlong,1.0D0)
 683        call reada(weightcard,'WSC',wsc,wlong)
 684        call reada(weightcard,'WSCP',wscp,wlong)
 685        call reada(weightcard,'WELEC',welec,1.0D0)
 686        call reada(weightcard,'WVDWPP',wvdwpp,welec)
 687        call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
 688        call reada(weightcard,'WCORR4',wcorr4,0.0D0)
 689        call reada(weightcard,'WCORR5',wcorr5,0.0D0)
 690        call reada(weightcard,'WCORR6',wcorr6,0.0D0)
 691        call reada(weightcard,'WTURN3',wturn3,1.0D0)
 692        call reada(weightcard,'WTURN4',wturn4,1.0D0)
 693        call reada(weightcard,'WTURN6',wturn6,1.0D0)
 694        call reada(weightcard,'WSCCOR',wsccor,1.0D0)
 695        call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
 696        call reada(weightcard,'WBOND',wbond,1.0D0)
 697        call reada(weightcard,'WTOR',wtor,1.0D0)
 698        call reada(weightcard,'WTORD',wtor_d,1.0D0)
 699        call reada(weightcard,'WANG',wang,1.0D0)
 700        call reada(weightcard,'WSCLOC',wscloc,1.0D0)
 701        call reada(weightcard,'SCAL14',scal14,0.4D0)
 702        call reada(weightcard,'SCALSCP',scalscp,1.0d0)
 703        call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
 704        call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
 705        call reada(weightcard,'TEMP0',temp0,300.0d0)
 706        call reada(weightcard,'WSHIELD',wshield,1.0d0)
 707        call reada(weightcard,'WLT',wliptran,0.0D0)
 708        call reada(weightcard,'WTUBE',wtube,1.0D0)
 709        if (index(weightcard,'SOFT').gt.0) ipot=6
 710 C 12/1/95 Added weight for the multi-body term WCORR
 711        call reada(weightcard,'WCORRH',wcorr,1.0D0)
 712        if (wcorr4.gt.0.0d0) wcorr=wcorr4
 713        weights(1)=wsc
 714        weights(2)=wscp
 715        weights(3)=welec
 716        weights(4)=wcorr
 717        weights(5)=wcorr5
 718        weights(6)=wcorr6
 719        weights(7)=wel_loc
 720        weights(8)=wturn3
 721        weights(9)=wturn4
 722        weights(10)=wturn6
 723        weights(11)=wang
 724        weights(12)=wscloc
 725        weights(13)=wtor
 726        weights(14)=wtor_d
 727        weights(15)=wstrain
 728        weights(16)=wvdwpp
 729        weights(17)=wbond
 730        weights(18)=scal14
 731        weights(21)=wsccor
 732       if(me.eq.king.or..not.out1file)
 733      & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
 734      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
 735      &  wturn4,wturn6
 736    10 format (/'Energy-term weights (unscaled):'//
 737      & 'WSCC=   ',f10.6,' (SC-SC)'/
 738      & 'WSCP=   ',f10.6,' (SC-p)'/
 739      & 'WELEC=  ',f10.6,' (p-p electr)'/
 740      & 'WVDWPP= ',f10.6,' (p-p VDW)'/
 741      & 'WBOND=  ',f10.6,' (stretching)'/
 742      & 'WANG=   ',f10.6,' (bending)'/
 743      & 'WSCLOC= ',f10.6,' (SC local)'/
 744      & 'WTOR=   ',f10.6,' (torsional)'/
 745      & 'WTORD=  ',f10.6,' (double torsional)'/
 746      & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
 747      & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
 748      & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
 749      & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
 750      & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
 751      & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
 752      & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
 753      & 'WTURN4= ',f10.6,' (turns, 4th order)'/
 754      & 'WTURN6= ',f10.6,' (turns, 6th order)')
 755       if(me.eq.king.or..not.out1file)then
 756        if (wcorr4.gt.0.0d0) then
 757         write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
 758      &   'between contact pairs of peptide groups'
 759         write (iout,'(2(a,f5.3/))')
 760      &  'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
 761      &  'Range of quenching the correlation terms:',2*delt_corr
 762        else if (wcorr.gt.0.0d0) then
 763         write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
 764      &   'between contact pairs of peptide groups'
 765        endif
 766        write (iout,'(a,f8.3)')
 767      &  'Scaling factor of 1,4 SC-p interactions:',scal14
 768        write (iout,'(a,f8.3)')
 769      &  'General scaling factor of SC-p interactions:',scalscp
 770       endif
 771       r0_corr=cutoff_corr-delt_corr
 772       do i=1,ntyp
 773         aad(i,1)=scalscp*aad(i,1)
 774         aad(i,2)=scalscp*aad(i,2)
 775         bad(i,1)=scalscp*bad(i,1)
 776         bad(i,2)=scalscp*bad(i,2)
 777       enddo
 778       call rescale_weights(t_bath)
 779       if(me.eq.king.or..not.out1file)
 780      & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
 781      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
 782      &  wturn4,wturn6
 783    22 format (/'Energy-term weights (scaled):'//
 784      & 'WSCC=   ',f10.6,' (SC-SC)'/
 785      & 'WSCP=   ',f10.6,' (SC-p)'/
 786      & 'WELEC=  ',f10.6,' (p-p electr)'/
 787      & 'WVDWPP= ',f10.6,' (p-p VDW)'/
 788      & 'WBOND=  ',f10.6,' (stretching)'/
 789      & 'WANG=   ',f10.6,' (bending)'/
 790      & 'WSCLOC= ',f10.6,' (SC local)'/
 791      & 'WTOR=   ',f10.6,' (torsional)'/
 792      & 'WTORD=  ',f10.6,' (double torsional)'/
 793      & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
 794      & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
 795      & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
 796      & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
 797      & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
 798      & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
 799      & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
 800      & 'WTURN4= ',f10.6,' (turns, 4th order)'/
 801      & 'WTURN6= ',f10.6,' (turns, 6th order)')
 802       if(me.eq.king.or..not.out1file)
 803      & write (iout,*) "Reference temperature for weights calculation:",
 804      &  temp0
 805       call reada(weightcard,"D0CM",d0cm,3.78d0)
 806       call reada(weightcard,"AKCM",akcm,15.1d0)
 807       call reada(weightcard,"AKTH",akth,11.0d0)
 808       call reada(weightcard,"AKCT",akct,12.0d0)
 809       call reada(weightcard,"V1SS",v1ss,-1.08d0)
 810       call reada(weightcard,"V2SS",v2ss,7.61d0)
 811       call reada(weightcard,"V3SS",v3ss,13.7d0)
 812       call reada(weightcard,"EBR",ebr,-5.50D0)
 813       call reada(weightcard,"ATRISS",atriss,0.301D0)
 814       call reada(weightcard,"BTRISS",btriss,0.021D0)
 815       call reada(weightcard,"CTRISS",ctriss,1.001D0)
 816       call reada(weightcard,"DTRISS",dtriss,1.001D0)
 817       call reada(weightcard,"LIPSCALE",lipscale,1.3D0)
 818       write (iout,*) "ATRISS=", atriss
 819       write (iout,*) "BTRISS=", btriss
 820       write (iout,*) "CTRISS=", ctriss
 821       write (iout,*) "DTRISS=", dtriss
 822        if (shield_mode.gt.0) then
 823         lipscale=0.0d0
 824         write (iout,*) "liscale not used in shield mode"
 825        endif
 826       write (iout,*) "lipid scaling factor", lipscale
 827       dyn_ss=(index(weightcard,'DYN_SS').gt.0)
 828       do i=1,maxres
 829         dyn_ss_mask(i)=.false.
 830       enddo
 831       do i=1,maxres-1
 832         do j=i+1,maxres
 833           dyn_ssbond_ij(i,j)=1.0d300
 834         enddo
 835       enddo
 836       call reada(weightcard,"HT",Ht,0.0D0)
 837       if (dyn_ss) then
 838         ss_depth=ebr/wsc-0.25*eps(1,1)
 839         Ht=Ht/wsc-0.25*eps(1,1)
 840         akcm=akcm/wsc
 841         akth=akth/wsc
 842         akct=akct/wsc
 843         v1ss=v1ss/wsc
 844         v2ss=v2ss/wsc
 845         v3ss=v3ss/wsc
 846       else
 847         if (wstrain.ne.0.0) then
 848          ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
 849         else
 850           ss_depth=0.0
 851         endif
 852       endif
 853       write (iout,*) "wshield,", wshield
 854       if(me.eq.king.or..not.out1file) then
 855        write (iout,*) "Parameters of the SS-bond potential:"
 856        write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
 857      & " AKCT",akct
 858        write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
 859        write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
 860        write (iout,*)" HT",Ht
 861        print *,'indpdb=',indpdb,' pdbref=',pdbref
 862       endif
 863       if (indpdb.gt.0 .or. pdbref) then
 864         read(inp,'(a)') pdbfile
 865         if(me.eq.king.or..not.out1file)
 866      &   write (iout,'(2a)') 'PDB data will be read from file ',
 867      &   pdbfile(:ilen(pdbfile))
 868         open(ipdbin,file=pdbfile,status='old',err=33)
 869         goto 34
 870   33    write (iout,'(a)') 'Error opening PDB file.'
 871         stop
 872   34    continue
 873 c        write (iout,*) 'Begin reading pdb data'
 874 c        call flush(iout)
 875         call readpdb
 876 c        write (iout,*) 'Finished reading pdb data'
 877 c        call flush(iout)
 878         if(me.eq.king.or..not.out1file)
 879      &   write (iout,'(a,i3,a,i3)')'nsup=',nsup,
 880      &   ' nstart_sup=',nstart_sup
 881         do i=1,nres
 882           itype_pdb(i)=itype(i)
 883         enddo
 884         close (ipdbin)
 885         nnt=nstart_sup
 886         nct=nstart_sup+nsup-1
 887         call contact(.false.,ncont_ref,icont_ref,co)
 888
 889         if (sideadd) then
 890          if(me.eq.king.or..not.out1file)
 891      &    write(iout,*)'Adding sidechains'
 892          maxsi=1000
 893          do i=2,nres-1
 894           iti=itype(i)
 895           if (iti.ne.10 .and. itype(i).ne.ntyp1) then
 896             nsi=0
 897             fail=.true.
 898             do while (fail.and.nsi.le.maxsi)
 899               call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
 900               nsi=nsi+1
 901             enddo
 902             if(fail) write(iout,*)'Adding sidechain failed for res ',
 903      &              i,' after ',nsi,' trials'
 904           endif
 905          enddo
 906         endif
 907       endif
 908       if (indpdb.eq.0) then
 909 C Read sequence if not taken from the pdb file.
 910         read (inp,*) nres
 911 c        print *,'nres=',nres
 912         if (iscode.gt.0) then
 913           read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
 914         else
 915           read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
 916         endif
 917 C Convert sequence to numeric code
 918         do i=1,nres
 919           itype(i)=rescode(i,sequence(i),iscode)
 920         enddo
 921 C Assign initial virtual bond lengths
 922         do i=2,nres
 923           vbld(i)=vbl
 924           vbld_inv(i)=vblinv
 925         enddo
 926         do i=2,nres-1
 927           vbld(i+nres)=dsc(iabs(itype(i)))
 928           vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
 929 c          write (iout,*) "i",i," itype",itype(i),
 930 c     &      " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
 931         enddo
 932       endif
 933 c      print *,nres
 934 c      print '(20i4)',(itype(i),i=1,nres)
 935       do i=1,nres
 936 #ifdef PROCOR
 937         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
 938 #else
 939         if (itype(i).eq.ntyp1) then
 940 #endif
 941           itel(i)=0
 942 #ifdef PROCOR
 943         else if (iabs(itype(i+1)).ne.20) then
 944 #else
 945         else if (iabs(itype(i)).ne.20) then
 946 #endif
 947           itel(i)=1
 948         else
 949           itel(i)=2
 950         endif
 951       enddo
 952       if(me.eq.king.or..not.out1file)then
 953        write (iout,*) "ITEL"
 954        do i=1,nres-1
 955          write (iout,*) i,itype(i),itel(i)
 956        enddo
 957        print *,'Call Read_Bridge.'
 958       endif
 959       call read_bridge
 960 C 8/13/98 Set limits to generating the dihedral angles
 961       do i=1,nres
 962         phibound(1,i)=-pi
 963         phibound(2,i)=pi
 964       enddo
 965       read (inp,*) ndih_constr
 966       if (ndih_constr.gt.0) then
 967 C        read (inp,*) ftors
 968         read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
 969      &  i=1,ndih_constr)
 970         if(me.eq.king.or..not.out1file)then
 971          write (iout,*)
 972      &   'There are',ndih_constr,' constraints on phi angles.'
 973          do i=1,ndih_constr
 974           write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
 975      &    ftors(i)
 976          enddo
 977         endif
 978         do i=1,ndih_constr
 979           phi0(i)=deg2rad*phi0(i)
 980           drange(i)=deg2rad*drange(i)
 981         enddo
 982 C        if(me.eq.king.or..not.out1file)
 983 C     &   write (iout,*) 'FTORS',ftors
 984         do i=1,ndih_constr
 985           ii = idih_constr(i)
 986           phibound(1,ii) = phi0(i)-drange(i)
 987           phibound(2,ii) = phi0(i)+drange(i)
 988         enddo
 989       endif
 990 C first setting the theta boundaries to 0 to pi
 991 C this mean that there is no energy penalty for any angle occuring this can be applied
 992 C for generate random conformation but is not implemented in this way
 993 C      do i=1,nres
 994 C        thetabound(1,i)=0
 995 C        thetabound(2,i)=pi
 996 C      enddo
 997 C begin reading theta constrains this is quartic constrains allowing to
 998 C have smooth second derivative
 999       if (with_theta_constr) then
1000 C with_theta_constr is keyword allowing for occurance of theta constrains
1001       read (inp,*) ntheta_constr
1002 C ntheta_constr is the number of theta constrains
1003       if (ntheta_constr.gt.0) then
1004 C        read (inp,*) ftors
1005         read (inp,*) (itheta_constr(i),theta_constr0(i),
1006      &  theta_drange(i),for_thet_constr(i),
1007      &  i=1,ntheta_constr)
1008 C the above code reads from 1 to ntheta_constr
1009 C itheta_constr(i) residue i for which is theta_constr
1010 C theta_constr0 the global minimum value
1011 C theta_drange is range for which there is no energy penalty
1012 C for_thet_constr is the force constant for quartic energy penalty
1013 C E=k*x**4
1014         if(me.eq.king.or..not.out1file)then
1015          write (iout,*)
1016      &   'There are',ntheta_constr,' constraints on phi angles.'
1017          do i=1,ntheta_constr
1018           write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
1019      &    theta_drange(i),
1020      &    for_thet_constr(i)
1021          enddo
1022         endif
1023         do i=1,ntheta_constr
1024           theta_constr0(i)=deg2rad*theta_constr0(i)
1025           theta_drange(i)=deg2rad*theta_drange(i)
1026         enddo
1027 C        if(me.eq.king.or..not.out1file)
1028 C     &   write (iout,*) 'FTORS',ftors
1029 C        do i=1,ntheta_constr
1030 C          ii = itheta_constr(i)
1031 C          thetabound(1,ii) = phi0(i)-drange(i)
1032 C          thetabound(2,ii) = phi0(i)+drange(i)
1033 C        enddo
1034       endif ! ntheta_constr.gt.0
1035       endif! with_theta_constr
1036 C
1037 C      with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
1038 C      write (iout,*) "with_dihed_constr ",with_dihed_constr
1039       nnt=1
1040 #ifdef MPI
1041       if (me.eq.king) then
1042 #endif
1043        write (iout,'(a)') 'Boundaries in phi angle sampling:'
1044        do i=1,nres
1045          write (iout,'(a3,i5,2f10.1)')
1046      &   restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
1047        enddo
1048 #ifdef MP
1049       endif
1050 #endif
1051       nct=nres
1052 cd      print *,'NNT=',NNT,' NCT=',NCT
1053       if (itype(1).eq.ntyp1) nnt=2
1054       if (itype(nres).eq.ntyp1) nct=nct-1
1055       if (pdbref) then
1056         if(me.eq.king.or..not.out1file)
1057      &   write (iout,'(a,i3)') 'nsup=',nsup
1058         nstart_seq=nnt
1059         if (nsup.le.(nct-nnt+1)) then
1060           do i=0,nct-nnt+1-nsup
1061             if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
1062               nstart_seq=nnt+i
1063               goto 111
1064             endif
1065           enddo
1066           write (iout,'(a)')
1067      &            'Error - sequences to be superposed do not match.'
1068           stop
1069         else
1070           do i=0,nsup-(nct-nnt+1)
1071             if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
1072      &      then
1073               nstart_sup=nstart_sup+i
1074               nsup=nct-nnt+1
1075               goto 111
1076             endif
1077           enddo
1078           write (iout,'(a)')
1079      &            'Error - sequences to be superposed do not match.'
1080         endif
1081   111   continue
1082         if (nsup.eq.0) nsup=nct-nnt
1083         if (nstart_sup.eq.0) nstart_sup=nnt
1084         if (nstart_seq.eq.0) nstart_seq=nnt
1085         if(me.eq.king.or..not.out1file)
1086      &   write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
1087      &                 ' nstart_seq=',nstart_seq
1088       endif
1089 c--- Zscore rms -------
1090       if (nz_start.eq.0) nz_start=nnt
1091       if (nz_end.eq.0 .and. nsup.gt.0) then
1092         nz_end=nnt+nsup-1
1093       else if (nz_end.eq.0) then
1094         nz_end=nct
1095       endif
1096       if(me.eq.king.or..not.out1file)then
1097        write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1098        write (iout,*) 'IZ_SC=',iz_sc
1099       endif
1100 c----------------------
1101       call init_int_table
1102       if (refstr) then
1103         if (.not.pdbref) then
1104           call read_angles(inp,*38)
1105           goto 39
1106    38     write (iout,'(a)') 'Error reading reference structure.'
1107 #ifdef MPI
1108           call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1109           stop 'Error reading reference structure'
1110 #endif
1111    39     call chainbuild_extconf
1112           call setup_var
1113 czscore          call geom_to_var(nvar,coord_exp_zs(1,1))
1114           nstart_sup=nnt
1115           nstart_seq=nnt
1116           nsup=nct-nnt+1
1117           kkk=1
1118           do i=1,2*nres
1119             do j=1,3
1120               cref(j,i,kkk)=c(j,i)
1121             enddo
1122           enddo
1123           call contact(.true.,ncont_ref,icont_ref,co)
1124         endif
1125         endif
1126         print *, "A TU"
1127         write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1128         call flush(iout)
1129         if (constr_dist.gt.0) call read_dist_constr
1130         write (iout,*) "After read_dist_constr nhpb",nhpb
1131         if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
1132         call hpb_partition
1133         if(me.eq.king.or..not.out1file)
1134      &   write (iout,*) 'Contact order:',co
1135         if (pdbref) then
1136         if(me.eq.king.or..not.out1file)
1137      &   write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1138         do i=1,ncont_ref
1139           do j=1,2
1140             icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1141           enddo
1142           if(me.eq.king.or..not.out1file)
1143      &     write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1144      &     icont_ref(1,i),' ',
1145      &     restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1146         enddo
1147         endif
1148 C      endif
1149       if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1150      &    .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1151      &    modecalc.ne.10) then
1152 C If input structure hasn't been supplied from the PDB file read or generate
1153 C initial geometry.
1154         if (iranconf.eq.0 .and. .not. extconf) then
1155           if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1156      &     write (iout,'(a)') 'Initial geometry will be read in.'
1157           if (read_cart) then
1158             read(inp,'(8f10.5)',end=36,err=36)
1159      &       ((c(l,k),l=1,3),k=1,nres),
1160      &       ((c(l,k+nres),l=1,3),k=nnt,nct)
1161             write (iout,*) "Exit READ_CART"
1162             write (iout,'(8f10.5)')
1163      &       ((c(l,k),l=1,3),k=1,nres),
1164      &       ((c(l,k+nres),l=1,3),k=nnt,nct)
1165             call int_from_cart1(.true.)
1166             write (iout,*) "Finish INT_TO_CART"
1167             do i=1,nres-1
1168               do j=1,3
1169                 dc(j,i)=c(j,i+1)-c(j,i)
1170                 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1171               enddo
1172             enddo
1173             do i=nnt,nct
1174               if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1175                 do j=1,3
1176                   dc(j,i+nres)=c(j,i+nres)-c(j,i)
1177                   dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1178                 enddo
1179               endif
1180             enddo
1181             return
1182           else
1183             call read_angles(inp,*36)
1184             call chainbuild_extconf
1185           endif
1186           goto 37
1187    36     write (iout,'(a)') 'Error reading angle file.'
1188 #ifdef MPI
1189           call mpi_finalize( MPI_COMM_WORLD,IERR )
1190 #endif
1191           stop 'Error reading angle file.'
1192    37     continue
1193         else if (extconf) then
1194          if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1195      &    write (iout,'(a)') 'Extended chain initial geometry.'
1196          do i=3,nres
1197           theta(i)=90d0*deg2rad
1198          enddo
1199          do i=4,nres
1200           phi(i)=180d0*deg2rad
1201          enddo
1202          do i=2,nres-1
1203           alph(i)=110d0*deg2rad
1204          enddo
1205          do i=2,nres-1
1206           omeg(i)=-120d0*deg2rad
1207           if (itype(i).le.0) omeg(i)=-omeg(i)
1208          enddo
1209         call chainbuild_extconf
1210         else
1211           if(me.eq.king.or..not.out1file)
1212      &     write (iout,'(a)') 'Random-generated initial geometry.'
1213
1214
1215 #ifdef MPI
1216           if (me.eq.king  .or. fg_rank.eq.0 .and. (
1217      &           modecalc.eq.12 .or. modecalc.eq.14) ) then
1218 #endif
1219             do itrial=1,100
1220               itmp=1
1221               call gen_rand_conf(itmp,*30)
1222               goto 40
1223    30         write (iout,*) 'Failed to generate random conformation',
1224      &          ', itrial=',itrial
1225               write (*,*) 'Processor:',me,
1226      &          ' Failed to generate random conformation',
1227      &          ' itrial=',itrial
1228               call intout
1229
1230 #ifdef AIX
1231               call flush_(iout)
1232 #else
1233               call flush(iout)
1234 #endif
1235             enddo
1236             write (iout,'(a,i3,a)') 'Processor:',me,
1237      &        ' error in generating random conformation.'
1238             write (*,'(a,i3,a)') 'Processor:',me,
1239      &        ' error in generating random conformation.'
1240             call flush(iout)
1241 #ifdef MPI
1242             call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1243    40       continue
1244           endif
1245 #else
1246           write (*,'(a)')
1247      &      ' error in generating random conformation.'
1248           stop
1249    40     continue
1250 #endif
1251         endif
1252       elseif (modecalc.eq.4) then
1253         read (inp,'(a)') intinname
1254         open (intin,file=intinname,status='old',err=333)
1255         if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1256      &  write (iout,'(a)') 'intinname',intinname
1257         write (*,'(a)') 'Processor',myrank,' intinname',intinname
1258         goto 334
1259   333   write (iout,'(2a)') 'Error opening angle file ',intinname
1260 #ifdef MPI
1261         call MPI_Finalize(MPI_COMM_WORLD,IERR)
1262 #endif
1263         stop 'Error opening angle file.'
1264   334   continue
1265
1266       endif
1267 C Generate distance constraints, if the PDB structure is to be regularized.
1268       if (nthread.gt.0) then
1269         call read_threadbase
1270       endif
1271       call setup_var
1272       if (me.eq.king .or. .not. out1file)
1273      & call intout
1274       if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1275         write (iout,'(/a,i3,a)')
1276      &  'The chain contains',ns,' disulfide-bridging cysteines.'
1277         write (iout,'(20i4)') (iss(i),i=1,ns)
1278        if (dyn_ss) then
1279           write(iout,*)"Running with dynamic disulfide-bond formation"
1280        else
1281         write (iout,'(/a/)') 'Pre-formed links are:'
1282         do i=1,nss
1283           i1=ihpb(i)-nres
1284           i2=jhpb(i)-nres
1285           it1=itype(i1)
1286           it2=itype(i2)
1287           write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1288      &    restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1289      &    ebr,forcon(i)
1290         enddo
1291         write (iout,'(a)')
1292        endif
1293       endif
1294       if (ns.gt.0.and.dyn_ss) then
1295           do i=nss+1,nhpb
1296             ihpb(i-nss)=ihpb(i)
1297             jhpb(i-nss)=jhpb(i)
1298             forcon(i-nss)=forcon(i)
1299             dhpb(i-nss)=dhpb(i)
1300           enddo
1301           nhpb=nhpb-nss
1302           nss=0
1303           call hpb_partition
1304           do i=1,ns
1305             dyn_ss_mask(iss(i))=.true.
1306           enddo
1307       endif
1308       if (i2ndstr.gt.0) call secstrp2dihc
1309 c      call geom_to_var(nvar,x)
1310 c      call etotal(energia(0))
1311 c      call enerprint(energia(0))
1312 c      call briefout(0,etot)
1313 c      stop
1314 cd    write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1315 cd    write (iout,'(a)') 'Variable list:'
1316 cd    write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1317 #ifdef MPI
1318       if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1319      &  write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1320      &  'Processor',myrank,': end reading molecular data.'
1321 #endif
1322       print *,"A TU?"
1323       return
1324       end
1325 c--------------------------------------------------------------------------
1326       logical function seq_comp(itypea,itypeb,length)
1327       implicit none
1328       integer length,itypea(length),itypeb(length)
1329       integer i
1330       do i=1,length
1331         if (itypea(i).ne.itypeb(i)) then
1332           seq_comp=.false.
1333           return
1334         endif
1335       enddo
1336       seq_comp=.true.
1337       return
1338       end
1339 c-----------------------------------------------------------------------------
1340       subroutine read_bridge
1341 C Read information about disulfide bridges.
1342       implicit real*8 (a-h,o-z)
1343       include 'DIMENSIONS'
1344 #ifdef MPI
1345       include 'mpif.h'
1346 #endif
1347       include 'COMMON.IOUNITS'
1348       include 'COMMON.GEO'
1349       include 'COMMON.VAR'
1350       include 'COMMON.INTERACT'
1351       include 'COMMON.LOCAL'
1352       include 'COMMON.NAMES'
1353       include 'COMMON.CHAIN'
1354       include 'COMMON.FFIELD'
1355       include 'COMMON.SBRIDGE'
1356       include 'COMMON.HEADER'
1357       include 'COMMON.CONTROL'
1358       include 'COMMON.DBASE'
1359       include 'COMMON.THREAD'
1360       include 'COMMON.TIME1'
1361       include 'COMMON.SETUP'
1362 C Read bridging residues.
1363       read (inp,*) ns,(iss(i),i=1,ns)
1364       print *,'ns=',ns
1365       if(me.eq.king.or..not.out1file)
1366      &  write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1367 C Check whether the specified bridging residues are cystines.
1368       do i=1,ns
1369         if (itype(iss(i)).ne.1) then
1370           if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1371      &   'Do you REALLY think that the residue ',
1372      &    restyp(itype(iss(i))),i,
1373      &   ' can form a disulfide bridge?!!!'
1374           write (*,'(2a,i3,a)')
1375      &   'Do you REALLY think that the residue ',
1376      &    restyp(itype(iss(i))),i,
1377      &   ' can form a disulfide bridge?!!!'
1378 #ifdef MPI
1379          call MPI_Finalize(MPI_COMM_WORLD,ierror)
1380          stop
1381 #endif
1382         endif
1383       enddo
1384 C Read preformed bridges.
1385       if (ns.gt.0) then
1386       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1387       if(fg_rank.eq.0)
1388      & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1389       if (nss.gt.0) then
1390         nhpb=nss
1391 C Check if the residues involved in bridges are in the specified list of
1392 C bridging residues.
1393         do i=1,nss
1394           do j=1,i-1
1395             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1396      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1397               write (iout,'(a,i3,a)') 'Disulfide pair',i,
1398      &      ' contains residues present in other pairs.'
1399               write (*,'(a,i3,a)') 'Disulfide pair',i,
1400      &      ' contains residues present in other pairs.'
1401 #ifdef MPI
1402               call MPI_Finalize(MPI_COMM_WORLD,ierror)
1403               stop
1404 #endif
1405             endif
1406           enddo
1407           do j=1,ns
1408             if (ihpb(i).eq.iss(j)) goto 10
1409           enddo
1410           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1411    10     continue
1412           do j=1,ns
1413             if (jhpb(i).eq.iss(j)) goto 20
1414           enddo
1415           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1416    20     continue
1417           dhpb(i)=dbr
1418           forcon(i)=fbr
1419         enddo
1420         do i=1,nss
1421           ihpb(i)=ihpb(i)+nres
1422           jhpb(i)=jhpb(i)+nres
1423         enddo
1424       endif
1425       endif
1426       return
1427       end
1428 c----------------------------------------------------------------------------
1429       subroutine read_x(kanal,*)
1430       implicit real*8 (a-h,o-z)
1431       include 'DIMENSIONS'
1432       include 'COMMON.GEO'
1433       include 'COMMON.VAR'
1434       include 'COMMON.CHAIN'
1435       include 'COMMON.IOUNITS'
1436       include 'COMMON.CONTROL'
1437       include 'COMMON.LOCAL'
1438       include 'COMMON.INTERACT'
1439 c Read coordinates from input
1440 c
1441       read(kanal,'(8f10.5)',end=10,err=10)
1442      &  ((c(l,k),l=1,3),k=1,nres),
1443      &  ((c(l,k+nres),l=1,3),k=nnt,nct)
1444       do j=1,3
1445         c(j,nres+1)=c(j,1)
1446         c(j,2*nres)=c(j,nres)
1447       enddo
1448       call int_from_cart1(.false.)
1449       do i=1,nres-1
1450         do j=1,3
1451           dc(j,i)=c(j,i+1)-c(j,i)
1452           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1453         enddo
1454       enddo
1455       do i=nnt,nct
1456         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1457           do j=1,3
1458             dc(j,i+nres)=c(j,i+nres)-c(j,i)
1459             dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1460           enddo
1461         endif
1462       enddo
1463
1464       return
1465    10 return1
1466       end
1467 c----------------------------------------------------------------------------
1468       subroutine read_threadbase
1469       implicit real*8 (a-h,o-z)
1470       include 'DIMENSIONS'
1471       include 'COMMON.IOUNITS'
1472       include 'COMMON.GEO'
1473       include 'COMMON.VAR'
1474       include 'COMMON.INTERACT'
1475       include 'COMMON.LOCAL'
1476       include 'COMMON.NAMES'
1477       include 'COMMON.CHAIN'
1478       include 'COMMON.FFIELD'
1479       include 'COMMON.SBRIDGE'
1480       include 'COMMON.HEADER'
1481       include 'COMMON.CONTROL'
1482       include 'COMMON.DBASE'
1483       include 'COMMON.THREAD'
1484       include 'COMMON.TIME1'
1485 C Read pattern database for threading.
1486       read (icbase,*) nseq
1487       do i=1,nseq
1488         read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1489      &   nres_base(2,i),nres_base(3,i)
1490         read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1491      &   nres_base(1,i))
1492 c       write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1493 c    &   nres_base(2,i),nres_base(3,i)
1494 c       write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1495 c    &   nres_base(1,i))
1496       enddo
1497       close (icbase)
1498       if (weidis.eq.0.0D0) weidis=0.1D0
1499       do i=nnt,nct
1500         do j=i+2,nct
1501           nhpb=nhpb+1
1502           ihpb(nhpb)=i
1503           jhpb(nhpb)=j
1504           forcon(nhpb)=weidis
1505         enddo
1506       enddo
1507       read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1508       write (iout,'(a,i5)') 'nexcl: ',nexcl
1509       write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1510       return
1511       end
1512 c------------------------------------------------------------------------------
1513       subroutine setup_var
1514       implicit real*8 (a-h,o-z)
1515       include 'DIMENSIONS'
1516       include 'COMMON.IOUNITS'
1517       include 'COMMON.GEO'
1518       include 'COMMON.VAR'
1519       include 'COMMON.INTERACT'
1520       include 'COMMON.LOCAL'
1521       include 'COMMON.NAMES'
1522       include 'COMMON.CHAIN'
1523       include 'COMMON.FFIELD'
1524       include 'COMMON.SBRIDGE'
1525       include 'COMMON.HEADER'
1526       include 'COMMON.CONTROL'
1527       include 'COMMON.DBASE'
1528       include 'COMMON.THREAD'
1529       include 'COMMON.TIME1'
1530 C Set up variable list.
1531       ntheta=nres-2
1532       nphi=nres-3
1533       nvar=ntheta+nphi
1534       nside=0
1535       do i=2,nres-1
1536         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1537           nside=nside+1
1538           ialph(i,1)=nvar+nside
1539           ialph(nside,2)=i
1540         endif
1541       enddo
1542       if (indphi.gt.0) then
1543         nvar=nphi
1544       else if (indback.gt.0) then
1545         nvar=nphi+ntheta
1546       else
1547         nvar=nvar+2*nside
1548       endif
1549 cd    write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1550       return
1551       end
1552 c----------------------------------------------------------------------------
1553       subroutine gen_dist_constr
1554 C Generate CA distance constraints.
1555       implicit real*8 (a-h,o-z)
1556       include 'DIMENSIONS'
1557       include 'COMMON.IOUNITS'
1558       include 'COMMON.GEO'
1559       include 'COMMON.VAR'
1560       include 'COMMON.INTERACT'
1561       include 'COMMON.LOCAL'
1562       include 'COMMON.NAMES'
1563       include 'COMMON.CHAIN'
1564       include 'COMMON.FFIELD'
1565       include 'COMMON.SBRIDGE'
1566       include 'COMMON.HEADER'
1567       include 'COMMON.CONTROL'
1568       include 'COMMON.DBASE'
1569       include 'COMMON.THREAD'
1570       include 'COMMON.TIME1'
1571       dimension itype_pdb(maxres)
1572       common /pizda/ itype_pdb
1573       character*2 iden
1574 cd      print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1575 cd      write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1576 cd     & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1577 cd     & ' nsup',nsup
1578       if (constr_dist.eq.11) then
1579              do i=nstart_sup,nstart_sup+nsup-1
1580               do j=i+2,nstart_sup+nsup-1
1581                  distance=dist(i,j)
1582                  if (distance.le.15.0) then
1583                  nhpb=nhpb+1
1584                  ihpb(nhpb)=i+nstart_seq-nstart_sup
1585                  jhpb(nhpb)=j+nstart_seq-nstart_sup
1586                  forcon(nhpb)=sqrt(0.04*distance)
1587                  fordepth(nhpb)=sqrt(40.0/distance)
1588                  dhpb(nhpb)=distance-0.1d0
1589                  dhpb1(nhpb)=distance+0.1d0
1590
1591 #ifdef MPI
1592           if (.not.out1file .or. me.eq.king)
1593      &    write (iout,'(a,3i5,f8.2,f10.2)') "+dist.constr ",
1594      &     nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1595 #else
1596           write (iout,'(a,3i5,f8.2,f10.2)') "+dist.constr ",
1597      &     nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1598 #endif
1599             endif
1600              enddo
1601            enddo
1602       else
1603       do i=nstart_sup,nstart_sup+nsup-1
1604 cd      write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1605 cd     &    ' seq_pdb', restyp(itype_pdb(i))
1606         do j=i+2,nstart_sup+nsup-1
1607           nhpb=nhpb+1
1608           ihpb(nhpb)=i+nstart_seq-nstart_sup
1609           jhpb(nhpb)=j+nstart_seq-nstart_sup
1610           forcon(nhpb)=weidis
1611           dhpb(nhpb)=dist(i,j)
1612         enddo
1613       enddo
1614       endif
1615 cd      write (iout,'(a)') 'Distance constraints:'
1616 cd      do i=nss+1,nhpb
1617 cd        ii=ihpb(i)
1618 cd        jj=jhpb(i)
1619 cd        iden='CA'
1620 cd        if (ii.gt.nres) then
1621 cd          iden='SC'
1622 cd          ii=ii-nres
1623 cd          jj=jj-nres
1624 cd        endif
1625 cd        write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1626 cd     &  restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1627 cd     &  dhpb(i),forcon(i)
1628 cd      enddo
1629       return
1630       end
1631 c----------------------------------------------------------------------------
1632       subroutine map_read
1633       implicit real*8 (a-h,o-z)
1634       include 'DIMENSIONS'
1635       include 'COMMON.MAP'
1636       include 'COMMON.IOUNITS'
1637       character*3 angid(4) /'THE','PHI','ALP','OME'/
1638       character*80 mapcard,ucase
1639       do imap=1,nmap
1640         read (inp,'(a)') mapcard
1641         mapcard=ucase(mapcard)
1642         if (index(mapcard,'PHI').gt.0) then
1643           kang(imap)=1
1644         else if (index(mapcard,'THE').gt.0) then
1645           kang(imap)=2
1646         else if (index(mapcard,'ALP').gt.0) then
1647           kang(imap)=3
1648         else if (index(mapcard,'OME').gt.0) then
1649           kang(imap)=4
1650         else
1651           write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1652           stop 'Error - illegal variable spec in MAP card.'
1653         endif
1654         call readi (mapcard,'RES1',res1(imap),0)
1655         call readi (mapcard,'RES2',res2(imap),0)
1656         if (res1(imap).eq.0) then
1657           res1(imap)=res2(imap)
1658         else if (res2(imap).eq.0) then
1659           res2(imap)=res1(imap)
1660         endif
1661         if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1662           write (iout,'(a)')
1663      &    'Error - illegal definition of variable group in MAP.'
1664           stop 'Error - illegal definition of variable group in MAP.'
1665         endif
1666         call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1667         call reada(mapcard,'TO',ang_to(imap),0.0D0)
1668         call readi(mapcard,'NSTEP',nstep(imap),0)
1669         if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1670           write (iout,'(a)')
1671      &     'Illegal boundary and/or step size specification in MAP.'
1672           stop 'Illegal boundary and/or step size specification in MAP.'
1673         endif
1674       enddo ! imap
1675       return
1676       end
1677 c----------------------------------------------------------------------------
1678       subroutine csaread
1679       implicit real*8 (a-h,o-z)
1680       include 'DIMENSIONS'
1681       include 'COMMON.IOUNITS'
1682       include 'COMMON.GEO'
1683       include 'COMMON.CSA'
1684       include 'COMMON.BANK'
1685       include 'COMMON.CONTROL'
1686       character*80 ucase
1687       character*620 mcmcard
1688       call card_concat(mcmcard)
1689
1690       call readi(mcmcard,'NCONF',nconf,50)
1691       call readi(mcmcard,'NADD',nadd,0)
1692       call readi(mcmcard,'JSTART',jstart,1)
1693       call readi(mcmcard,'JEND',jend,1)
1694       call readi(mcmcard,'NSTMAX',nstmax,500000)
1695       call readi(mcmcard,'N0',n0,1)
1696       call readi(mcmcard,'N1',n1,6)
1697       call readi(mcmcard,'N2',n2,4)
1698       call readi(mcmcard,'N3',n3,0)
1699       call readi(mcmcard,'N4',n4,0)
1700       call readi(mcmcard,'N5',n5,0)
1701       call readi(mcmcard,'N6',n6,10)
1702       call readi(mcmcard,'N7',n7,0)
1703       call readi(mcmcard,'N8',n8,0)
1704       call readi(mcmcard,'N9',n9,0)
1705       call readi(mcmcard,'N14',n14,0)
1706       call readi(mcmcard,'N15',n15,0)
1707       call readi(mcmcard,'N16',n16,0)
1708       call readi(mcmcard,'N17',n17,0)
1709       call readi(mcmcard,'N18',n18,0)
1710
1711       vdisulf=(index(mcmcard,'DYNSS').gt.0)
1712
1713       call readi(mcmcard,'NDIFF',ndiff,2)
1714       call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1715       call readi(mcmcard,'IS1',is1,1)
1716       call readi(mcmcard,'IS2',is2,8)
1717       call readi(mcmcard,'NRAN0',nran0,4)
1718       call readi(mcmcard,'NRAN1',nran1,2)
1719       call readi(mcmcard,'IRR',irr,1)
1720       call readi(mcmcard,'NSEED',nseed,20)
1721       call readi(mcmcard,'NTOTAL',ntotal,10000)
1722       call reada(mcmcard,'CUT1',cut1,2.0d0)
1723       call reada(mcmcard,'CUT2',cut2,5.0d0)
1724       call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1725       call readi(mcmcard,'ICMAX',icmax,3)
1726       call readi(mcmcard,'IRESTART',irestart,0)
1727 c!bankt      call readi(mcmcard,'NBANKTM',ntbankm,0)
1728       ntbankm=0
1729 c!bankt
1730       call reada(mcmcard,'DELE',dele,20.0d0)
1731       call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1732       call readi(mcmcard,'IREF',iref,0)
1733       call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1734       call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1735       call readi(mcmcard,'NCONF_IN',nconf_in,0)
1736       call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1737       write (iout,*) "NCONF_IN",nconf_in
1738       return
1739       end
1740 c----------------------------------------------------------------------------
1741 cfmc      subroutine mcmfread
1742 cfmc      implicit real*8 (a-h,o-z)
1743 cfmc      include 'DIMENSIONS'
1744 cfmc      include 'COMMON.MCMF'
1745 cfmc      include 'COMMON.IOUNITS'
1746 cfmc      include 'COMMON.GEO'
1747 cfmc      character*80 ucase
1748 cfmc      character*620 mcmcard
1749 cfmc      call card_concat(mcmcard)
1750 cfmc
1751 cfmc      call readi(mcmcard,'MAXRANT',maxrant,1000)
1752 cfmc      write(iout,*)'MAXRANT=',maxrant
1753 cfmc      call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1754 cfmc      write(iout,*)'MAXFAM=',maxfam
1755 cfmc      call readi(mcmcard,'NNET1',nnet1,5)
1756 cfmc      write(iout,*)'NNET1=',nnet1
1757 cfmc      call readi(mcmcard,'NNET2',nnet2,4)
1758 cfmc      write(iout,*)'NNET2=',nnet2
1759 cfmc      call readi(mcmcard,'NNET3',nnet3,4)
1760 cfmc      write(iout,*)'NNET3=',nnet3
1761 cfmc      call readi(mcmcard,'ILASTT',ilastt,0)
1762 cfmc      write(iout,*)'ILASTT=',ilastt
1763 cfmc      call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1764 cfmc      write(iout,*)'MAXSTR=',maxstr
1765 cfmc      maxstr_f=maxstr/maxfam
1766 cfmc      write(iout,*)'MAXSTR_F=',maxstr_f
1767 cfmc      call readi(mcmcard,'NMCMF',nmcmf,10)
1768 cfmc      write(iout,*)'NMCMF=',nmcmf
1769 cfmc      call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1770 cfmc      write(iout,*)'IFOCUS=',ifocus
1771 cfmc      call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1772 cfmc      write(iout,*)'NLOCMCMF=',nlocmcmf
1773 cfmc      call readi(mcmcard,'INTPRT',intprt,1000)
1774 cfmc      write(iout,*)'INTPRT=',intprt
1775 cfmc      call readi(mcmcard,'IPRT',iprt,100)
1776 cfmc      write(iout,*)'IPRT=',iprt
1777 cfmc      call readi(mcmcard,'IMAXTR',imaxtr,100)
1778 cfmc      write(iout,*)'IMAXTR=',imaxtr
1779 cfmc      call readi(mcmcard,'MAXEVEN',maxeven,1000)
1780 cfmc      write(iout,*)'MAXEVEN=',maxeven
1781 cfmc      call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1782 cfmc      write(iout,*)'MAXEVEN1=',maxeven1
1783 cfmc      call readi(mcmcard,'INIMIN',inimin,200)
1784 cfmc      write(iout,*)'INIMIN=',inimin
1785 cfmc      call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1786 cfmc      write(iout,*)'NSTEPMCMF=',nstepmcmf
1787 cfmc      call readi(mcmcard,'NTHREAD',nthread,5)
1788 cfmc      write(iout,*)'NTHREAD=',nthread
1789 cfmc      call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1790 cfmc      write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1791 cfmc      call readi(mcmcard,'MAXPERT',maxpert,9)
1792 cfmc      write(iout,*)'MAXPERT=',maxpert
1793 cfmc      call readi(mcmcard,'IRMSD',irmsd,1)
1794 cfmc      write(iout,*)'IRMSD=',irmsd
1795 cfmc      call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1796 cfmc      write(iout,*)'DENEMIN=',denemin
1797 cfmc      call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1798 cfmc      write(iout,*)'RCUT1S=',rcut1s
1799 cfmc      call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1800 cfmc      write(iout,*)'RCUT1E=',rcut1e
1801 cfmc      call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1802 cfmc      write(iout,*)'RCUT2S=',rcut2s
1803 cfmc      call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1804 cfmc      write(iout,*)'RCUT2E=',rcut2e
1805 cfmc      call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1806 cfmc      write(iout,*)'DPERT1=',d_pert1
1807 cfmc      call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1808 cfmc      write(iout,*)'DPERT1A=',d_pert1a
1809 cfmc      call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1810 cfmc      write(iout,*)'DPERT2=',d_pert2
1811 cfmc      call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1812 cfmc      write(iout,*)'DPERT2A=',d_pert2a
1813 cfmc      call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1814 cfmc      write(iout,*)'DPERT2B=',d_pert2b
1815 cfmc      call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1816 cfmc      write(iout,*)'DPERT2C=',d_pert2c
1817 cfmc      d_pert1=deg2rad*d_pert1
1818 cfmc      d_pert1a=deg2rad*d_pert1a
1819 cfmc      d_pert2=deg2rad*d_pert2
1820 cfmc      d_pert2a=deg2rad*d_pert2a
1821 cfmc      d_pert2b=deg2rad*d_pert2b
1822 cfmc      d_pert2c=deg2rad*d_pert2c
1823 cfmc      call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1824 cfmc      write(iout,*)'KT_MCMF1=',kt_mcmf1
1825 cfmc      call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1826 cfmc      write(iout,*)'KT_MCMF2=',kt_mcmf2
1827 cfmc      call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1828 cfmc      write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1829 cfmc      call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1830 cfmc      write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1831 cfmc      call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1832 cfmc      write(iout,*)'RCUTINI=',rcutini
1833 cfmc      call reada(mcmcard,'GRAT',grat,0.5D0)
1834 cfmc      write(iout,*)'GRAT=',grat
1835 cfmc      call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1836 cfmc      write(iout,*)'BIAS_MCMF=',bias_mcmf
1837 cfmc
1838 cfmc      return
1839 cfmc      end
1840 c----------------------------------------------------------------------------
1841       subroutine mcmread
1842       implicit real*8 (a-h,o-z)
1843       include 'DIMENSIONS'
1844       include 'COMMON.MCM'
1845       include 'COMMON.MCE'
1846       include 'COMMON.IOUNITS'
1847       character*80 ucase
1848       character*320 mcmcard
1849       call card_concat(mcmcard)
1850       call readi(mcmcard,'MAXACC',maxacc,100)
1851       call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1852       call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1853       call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1854       call readi(mcmcard,'MAXREPM',maxrepm,200)
1855       call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1856       call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1857       call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1858       call reada(mcmcard,'E_UP',e_up,5.0D0)
1859       call reada(mcmcard,'DELTE',delte,0.1D0)
1860       call readi(mcmcard,'NSWEEP',nsweep,5)
1861       call readi(mcmcard,'NSTEPH',nsteph,0)
1862       call readi(mcmcard,'NSTEPC',nstepc,0)
1863       call reada(mcmcard,'TMIN',tmin,298.0D0)
1864       call reada(mcmcard,'TMAX',tmax,298.0D0)
1865       call readi(mcmcard,'NWINDOW',nwindow,0)
1866       call readi(mcmcard,'PRINT_MC',print_mc,0)
1867       print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1868       print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1869       ent_read=(index(mcmcard,'ENT_READ').gt.0)
1870       call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1871       call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1872       call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1873       call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1874       call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1875       if (nwindow.gt.0) then
1876         read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1877         do i=1,nwindow
1878           winlen(i)=winend(i)-winstart(i)+1
1879         enddo
1880       endif
1881       if (tmax.lt.tmin) tmax=tmin
1882       if (tmax.eq.tmin) then
1883         nstepc=0
1884         nsteph=0
1885       endif
1886       if (nstepc.gt.0 .and. nsteph.gt.0) then
1887         tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1888         tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1889       endif
1890 C Probabilities of different move types
1891       sumpro_type(0)=0.0D0
1892       call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1893       call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1894       sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1895       call reada(mcmcard,'THETA'     ,sumpro_type(3),0.0d0)
1896       sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1897       call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1898       sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1899       do i=1,MaxMoveType
1900         print *,'i',i,' sumprotype',sumpro_type(i)
1901         sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1902         print *,'i',i,' sumprotype',sumpro_type(i)
1903       enddo
1904       return
1905       end
1906 c----------------------------------------------------------------------------
1907       subroutine read_minim
1908       implicit real*8 (a-h,o-z)
1909       include 'DIMENSIONS'
1910       include 'COMMON.MINIM'
1911       include 'COMMON.IOUNITS'
1912       character*80 ucase
1913       character*320 minimcard
1914       call card_concat(minimcard)
1915       call readi(minimcard,'MAXMIN',maxmin,2000)
1916       call readi(minimcard,'MAXFUN',maxfun,5000)
1917       call readi(minimcard,'MINMIN',minmin,maxmin)
1918       call readi(minimcard,'MINFUN',minfun,maxmin)
1919       call reada(minimcard,'TOLF',tolf,1.0D-2)
1920       call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1921       print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1922       print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1923       print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1924       write (iout,'(/80(1h*)/20x,a/80(1h*))')
1925      &         'Options in energy minimization:'
1926       write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1927      & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1928      & 'MinMin:',MinMin,' MinFun:',MinFun,
1929      & ' TolF:',TolF,' RTolF:',RTolF
1930       return
1931       end
1932 c----------------------------------------------------------------------------
1933       subroutine read_angles(kanal,*)
1934       implicit real*8 (a-h,o-z)
1935       include 'DIMENSIONS'
1936       include 'COMMON.GEO'
1937       include 'COMMON.VAR'
1938       include 'COMMON.CHAIN'
1939       include 'COMMON.IOUNITS'
1940       include 'COMMON.CONTROL'
1941 c Read angles from input
1942 c
1943        read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1944        read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1945        read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1946        read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1947
1948        do i=1,nres
1949 c 9/7/01 avoid 180 deg valence angle
1950         if (theta(i).gt.179.99d0) theta(i)=179.99d0
1951 c
1952         theta(i)=deg2rad*theta(i)
1953         phi(i)=deg2rad*phi(i)
1954         alph(i)=deg2rad*alph(i)
1955         omeg(i)=deg2rad*omeg(i)
1956        enddo
1957       return
1958    10 return1
1959       end
1960 c----------------------------------------------------------------------------
1961       subroutine reada(rekord,lancuch,wartosc,default)
1962       implicit none
1963       character*(*) rekord,lancuch
1964       double precision wartosc,default
1965       integer ilen,iread
1966       external ilen
1967       iread=index(rekord,lancuch)
1968       if (iread.eq.0) then
1969         wartosc=default
1970         return
1971       endif
1972       iread=iread+ilen(lancuch)+1
1973       read (rekord(iread:),*,err=10,end=10) wartosc
1974       return
1975   10  wartosc=default
1976       return
1977       end
1978 c----------------------------------------------------------------------------
1979       subroutine readi(rekord,lancuch,wartosc,default)
1980       implicit none
1981       character*(*) rekord,lancuch
1982       integer wartosc,default
1983       integer ilen,iread
1984       external ilen
1985       iread=index(rekord,lancuch)
1986       if (iread.eq.0) then
1987         wartosc=default
1988         return
1989       endif
1990       iread=iread+ilen(lancuch)+1
1991       read (rekord(iread:),*,err=10,end=10) wartosc
1992       return
1993   10  wartosc=default
1994       return
1995       end
1996 c----------------------------------------------------------------------------
1997       subroutine multreadi(rekord,lancuch,tablica,dim,default)
1998       implicit none
1999       integer dim,i
2000       integer tablica(dim),default
2001       character*(*) rekord,lancuch
2002       character*80 aux
2003       integer ilen,iread
2004       external ilen
2005       do i=1,dim
2006         tablica(i)=default
2007       enddo
2008       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2009       if (iread.eq.0) return
2010       iread=iread+ilen(lancuch)+1
2011
2012       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2013    10 return
2014       end
2015 c----------------------------------------------------------------------------
2016       subroutine multreada(rekord,lancuch,tablica,dim,default)
2017       implicit none
2018       integer dim,i
2019       double precision tablica(dim),default
2020       character*(*) rekord,lancuch
2021       character*80 aux
2022       integer ilen,iread
2023       external ilen
2024       do i=1,dim
2025         tablica(i)=default
2026       enddo
2027       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2028       if (iread.eq.0) return
2029       iread=iread+ilen(lancuch)+1
2030       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2031    10 return
2032       end
2033 c----------------------------------------------------------------------------
2034       subroutine openunits
2035       implicit real*8 (a-h,o-z)
2036       include 'DIMENSIONS'
2037 #ifdef MPI
2038       include 'mpif.h'
2039       character*16 form,nodename
2040       integer nodelen
2041 #endif
2042       include 'COMMON.SETUP'
2043       include 'COMMON.IOUNITS'
2044       include 'COMMON.MD'
2045       include 'COMMON.CONTROL'
2046       integer lenpre,lenpot,ilen,lentmp
2047       external ilen
2048       character*3 out1file_text,ucase
2049       character*3 ll
2050       external ucase
2051 c      print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
2052       call getenv_loc("PREFIX",prefix)
2053       pref_orig = prefix
2054       call getenv_loc("POT",pot)
2055       call getenv_loc("DIRTMP",tmpdir)
2056       call getenv_loc("CURDIR",curdir)
2057       call getenv_loc("OUT1FILE",out1file_text)
2058 c      print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
2059       out1file_text=ucase(out1file_text)
2060       if (out1file_text(1:1).eq."Y") then
2061         out1file=.true.
2062       else
2063         out1file=fg_rank.gt.0
2064       endif
2065       lenpre=ilen(prefix)
2066       lenpot=ilen(pot)
2067       lentmp=ilen(tmpdir)
2068       if (lentmp.gt.0) then
2069           write (*,'(80(1h!))')
2070           write (*,'(a,19x,a,19x,a)') "!","  A T T E N T I O N  ","!"
2071           write (*,'(80(1h!))')
2072           write (*,*)"All output files will be on node /tmp directory."
2073 #ifdef MPI
2074         call  MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
2075         if (me.eq.king) then
2076           write (*,*) "The master node is ",nodename
2077         else if (fg_rank.eq.0) then
2078           write (*,*) "I am the CG slave node ",nodename
2079         else
2080           write (*,*) "I am the FG slave node ",nodename
2081         endif
2082 #endif
2083         PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
2084         lenpre = lentmp+lenpre+1
2085       endif
2086       entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
2087 C Get the names and open the input files
2088 #if defined(WINIFL) || defined(WINPGI)
2089       open(1,file=pref_orig(:ilen(pref_orig))//
2090      &  '.inp',status='old',readonly,shared)
2091        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2092 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2093 C Get parameter filenames and open the parameter files.
2094       call getenv_loc('BONDPAR',bondname)
2095       open (ibond,file=bondname,status='old',readonly,shared)
2096       call getenv_loc('THETPAR',thetname)
2097       open (ithep,file=thetname,status='old',readonly,shared)
2098       call getenv_loc('ROTPAR',rotname)
2099       open (irotam,file=rotname,status='old',readonly,shared)
2100       call getenv_loc('TORPAR',torname)
2101       open (itorp,file=torname,status='old',readonly,shared)
2102       call getenv_loc('TORDPAR',tordname)
2103       open (itordp,file=tordname,status='old',readonly,shared)
2104       call getenv_loc('TORKCC',torkccname)
2105       open (itorkcc,file=torkccname,status='old',readonly,shared)
2106       call getenv_loc('THETKCC',thetkccname)
2107       open (ithetkcc,file=thetkccname,status='old',readonly,shared)
2108       call getenv_loc('FOURIER',fouriername)
2109       open (ifourier,file=fouriername,status='old',readonly,shared)
2110       call getenv_loc('ELEPAR',elename)
2111       open (ielep,file=elename,status='old',readonly,shared)
2112       call getenv_loc('SIDEPAR',sidename)
2113       open (isidep,file=sidename,status='old',readonly,shared)
2114       call getenv_loc('LIPTRANPAR',liptranname)
2115       open (iliptranpar,file=liptranname,status='old',readonly,shared)
2116 #elif (defined CRAY) || (defined AIX)
2117       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2118      &  action='read')
2119 c      print *,"Processor",myrank," opened file 1"
2120       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2121 c      print *,"Processor",myrank," opened file 9"
2122 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2123 C Get parameter filenames and open the parameter files.
2124       call getenv_loc('BONDPAR',bondname)
2125       open (ibond,file=bondname,status='old',action='read')
2126 c      print *,"Processor",myrank," opened file IBOND"
2127       call getenv_loc('THETPAR',thetname)
2128       open (ithep,file=thetname,status='old',action='read')
2129 c      print *,"Processor",myrank," opened file ITHEP"
2130       call getenv_loc('ROTPAR',rotname)
2131       open (irotam,file=rotname,status='old',action='read')
2132 c      print *,"Processor",myrank," opened file IROTAM"
2133       call getenv_loc('TORPAR',torname)
2134       open (itorp,file=torname,status='old',action='read')
2135 c      print *,"Processor",myrank," opened file ITORP"
2136       call getenv_loc('TORDPAR',tordname)
2137       open (itordp,file=tordname,status='old',action='read')
2138       call getenv_loc('TORKCC',torkccname)
2139       open (itorkcc,file=torkccname,status='old',action='read')
2140       call getenv_loc('THETKCC',thetkccname)
2141       open (ithetkcc,file=thetkccname,status='old',action='read')
2142 c      print *,"Processor",myrank," opened file ITORDP"
2143       call getenv_loc('SCCORPAR',sccorname)
2144       open (isccor,file=sccorname,status='old',action='read')
2145 c      print *,"Processor",myrank," opened file ISCCOR"
2146       call getenv_loc('FOURIER',fouriername)
2147       open (ifourier,file=fouriername,status='old',action='read')
2148 c      print *,"Processor",myrank," opened file IFOURIER"
2149       call getenv_loc('ELEPAR',elename)
2150       open (ielep,file=elename,status='old',action='read')
2151 c      print *,"Processor",myrank," opened file IELEP"
2152       call getenv_loc('SIDEPAR',sidename)
2153       open (isidep,file=sidename,status='old',action='read')
2154       call getenv_loc('LIPTRANPAR',liptranname)
2155       open (iliptranpar,file=liptranname,status='old',action='read')
2156 c      print *,"Processor",myrank," opened file ISIDEP"
2157 c      print *,"Processor",myrank," opened parameter files"
2158 #elif (defined G77)
2159       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2160       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2161 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2162 C Get parameter filenames and open the parameter files.
2163       call getenv_loc('BONDPAR',bondname)
2164       open (ibond,file=bondname,status='old')
2165       call getenv_loc('THETPAR',thetname)
2166       open (ithep,file=thetname,status='old')
2167       call getenv_loc('ROTPAR',rotname)
2168       open (irotam,file=rotname,status='old')
2169       call getenv_loc('TORPAR',torname)
2170       open (itorp,file=torname,status='old')
2171       call getenv_loc('TORDPAR',tordname)
2172       open (itordp,file=tordname,status='old')
2173       call getenv_loc('TORKCC',torkccname)
2174       open (itorkcc,file=torkccname,status='old')
2175       call getenv_loc('THETKCC',thetkccname)
2176       open (ithetkcc,file=thetkccname,status='old')
2177       call getenv_loc('SCCORPAR',sccorname)
2178       open (isccor,file=sccorname,status='old')
2179       call getenv_loc('FOURIER',fouriername)
2180       open (ifourier,file=fouriername,status='old')
2181       call getenv_loc('ELEPAR',elename)
2182       open (ielep,file=elename,status='old')
2183       call getenv_loc('SIDEPAR',sidename)
2184       open (isidep,file=sidename,status='old')
2185       call getenv_loc('LIPTRANPAR',liptranname)
2186       open (iliptranpar,file=liptranname,status='old')
2187 #else
2188       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2189      &  readonly)
2190        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2191 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2192 C Get parameter filenames and open the parameter files.
2193       call getenv_loc('BONDPAR',bondname)
2194       open (ibond,file=bondname,status='old',readonly)
2195       call getenv_loc('THETPAR',thetname)
2196       open (ithep,file=thetname,status='old',readonly)
2197       call getenv_loc('ROTPAR',rotname)
2198       open (irotam,file=rotname,status='old',readonly)
2199       call getenv_loc('TORPAR',torname)
2200       open (itorp,file=torname,status='old',readonly)
2201       call getenv_loc('TORDPAR',tordname)
2202       open (itordp,file=tordname,status='old',readonly)
2203       call getenv_loc('TORKCC',torkccname)
2204       open (itorkcc,file=torkccname,status='old',readonly)
2205       call getenv_loc('THETKCC',thetkccname)
2206       open (ithetkcc,file=thetkccname,status='old',readonly)
2207       call getenv_loc('SCCORPAR',sccorname)
2208       open (isccor,file=sccorname,status='old',readonly)
2209 #ifndef CRYST_THETA
2210       call getenv_loc('THETPARPDB',thetname_pdb)
2211       print *,"thetname_pdb ",thetname_pdb
2212       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2213       print *,ithep_pdb," opened"
2214 #endif
2215       call getenv_loc('FOURIER',fouriername)
2216       open (ifourier,file=fouriername,status='old',readonly)
2217       call getenv_loc('ELEPAR',elename)
2218       open (ielep,file=elename,status='old',readonly)
2219       call getenv_loc('SIDEPAR',sidename)
2220       open (isidep,file=sidename,status='old',readonly)
2221       call getenv_loc('LIPTRANPAR',liptranname)
2222       open (iliptranpar,file=liptranname,status='old',action='read')
2223 #ifndef CRYST_SC
2224       call getenv_loc('ROTPARPDB',rotname_pdb)
2225       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2226 #endif
2227 #endif
2228       call getenv_loc('TUBEPAR',tubename)
2229       open (itube,file=tubename,status='old',readonly)
2230 #ifndef OLDSCP
2231 C
2232 C 8/9/01 In the newest version SCp interaction constants are read from a file
2233 C Use -DOLDSCP to use hard-coded constants instead.
2234 C
2235       call getenv_loc('SCPPAR',scpname)
2236 #if defined(WINIFL) || defined(WINPGI)
2237       open (iscpp,file=scpname,status='old',readonly,shared)
2238 #elif (defined CRAY)  || (defined AIX)
2239       open (iscpp,file=scpname,status='old',action='read')
2240 #elif (defined G77)
2241       open (iscpp,file=scpname,status='old')
2242 #else
2243       open (iscpp,file=scpname,status='old',readonly)
2244 #endif
2245 #endif
2246       call getenv_loc('PATTERN',patname)
2247 #if defined(WINIFL) || defined(WINPGI)
2248       open (icbase,file=patname,status='old',readonly,shared)
2249 #elif (defined CRAY)  || (defined AIX)
2250       open (icbase,file=patname,status='old',action='read')
2251 #elif (defined G77)
2252       open (icbase,file=patname,status='old')
2253 #else
2254       open (icbase,file=patname,status='old',readonly)
2255 #endif
2256 #ifdef MPI
2257 C Open output file only for CG processes
2258 c      print *,"Processor",myrank," fg_rank",fg_rank
2259       if (fg_rank.eq.0) then
2260
2261       if (nodes.eq.1) then
2262         npos=3
2263       else
2264         npos = dlog10(dfloat(nodes-1))+1
2265       endif
2266       if (npos.lt.3) npos=3
2267       write (liczba,'(i1)') npos
2268       form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2269      &  //')'
2270       write (liczba,form) me
2271       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2272      &  liczba(:ilen(liczba))
2273       intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2274      &  //'.int'
2275       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2276      &  //'.pdb'
2277       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2278      &  liczba(:ilen(liczba))//'.mol2'
2279       statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2280      &  liczba(:ilen(liczba))//'.stat'
2281       if (lentmp.gt.0)
2282      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2283      &      //liczba(:ilen(liczba))//'.stat')
2284       rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2285      &  //'.rst'
2286       if(usampl) then
2287           qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2288      & liczba(:ilen(liczba))//'.const'
2289       endif
2290
2291       endif
2292 #else
2293       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2294       intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2295       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2296       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2297       statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2298       if (lentmp.gt.0)
2299      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2300      &    //'.stat')
2301       rest2name=prefix(:ilen(prefix))//'.rst'
2302       if(usampl) then
2303          qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2304       endif
2305 #endif
2306 #if defined(AIX) || defined(PGI)
2307       if (me.eq.king .or. .not. out1file)
2308      &   open(iout,file=outname,status='unknown')
2309 #ifdef DEBUG
2310       if (fg_rank.gt.0) then
2311         write (liczba,'(i3.3)') myrank/nfgtasks
2312         write (ll,'(bz,i3.3)') fg_rank
2313         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2314      &   status='unknown')
2315       endif
2316 #endif
2317       if(me.eq.king) then
2318        open(igeom,file=intname,status='unknown',position='append')
2319        open(ipdb,file=pdbname,status='unknown')
2320        open(imol2,file=mol2name,status='unknown')
2321        open(istat,file=statname,status='unknown',position='append')
2322       else
2323 c1out       open(iout,file=outname,status='unknown')
2324       endif
2325 #else
2326       if (me.eq.king .or. .not.out1file)
2327      &    open(iout,file=outname,status='unknown')
2328 #ifdef DEBUG
2329       if (fg_rank.gt.0) then
2330         write (liczba,'(i3.3)') myrank/nfgtasks
2331         write (ll,'(bz,i3.3)') fg_rank
2332         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2333      &   status='unknown')
2334       endif
2335 #endif
2336       if(me.eq.king) then
2337        open(igeom,file=intname,status='unknown',access='append')
2338        open(ipdb,file=pdbname,status='unknown')
2339        open(imol2,file=mol2name,status='unknown')
2340        open(istat,file=statname,status='unknown',access='append')
2341       else
2342 c1out       open(iout,file=outname,status='unknown')
2343       endif
2344 #endif
2345       csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2346       csa_seed=prefix(:lenpre)//'.CSA.seed'
2347       csa_history=prefix(:lenpre)//'.CSA.history'
2348       csa_bank=prefix(:lenpre)//'.CSA.bank'
2349       csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2350       csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2351       csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2352 c!bankt      csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2353       csa_int=prefix(:lenpre)//'.int'
2354       csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2355       csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2356       csa_in=prefix(:lenpre)//'.CSA.in'
2357 c      print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2358 C Write file names
2359       if (me.eq.king)then
2360       write (iout,'(80(1h-))')
2361       write (iout,'(30x,a)') "FILE ASSIGNMENT"
2362       write (iout,'(80(1h-))')
2363       write (iout,*) "Input file                      : ",
2364      &  pref_orig(:ilen(pref_orig))//'.inp'
2365       write (iout,*) "Output file                     : ",
2366      &  outname(:ilen(outname))
2367       write (iout,*)
2368       write (iout,*) "Sidechain potential file        : ",
2369      &  sidename(:ilen(sidename))
2370 #ifndef OLDSCP
2371       write (iout,*) "SCp potential file              : ",
2372      &  scpname(:ilen(scpname))
2373 #endif
2374       write (iout,*) "Electrostatic potential file    : ",
2375      &  elename(:ilen(elename))
2376       write (iout,*) "Cumulant coefficient file       : ",
2377      &  fouriername(:ilen(fouriername))
2378       write (iout,*) "Torsional parameter file        : ",
2379      &  torname(:ilen(torname))
2380       write (iout,*) "Double torsional parameter file : ",
2381      &  tordname(:ilen(tordname))
2382       write (iout,*) "SCCOR parameter file : ",
2383      &  sccorname(:ilen(sccorname))
2384       write (iout,*) "Bond & inertia constant file    : ",
2385      &  bondname(:ilen(bondname))
2386       write (iout,*) "Bending parameter file          : ",
2387      &  thetname(:ilen(thetname))
2388       write (iout,*) "Rotamer parameter file          : ",
2389      &  rotname(:ilen(rotname))
2390       write (iout,*) "Thetpdb parameter file          : ",
2391      &  thetname_pdb(:ilen(thetname_pdb))
2392       write (iout,*) "Threading database              : ",
2393      &  patname(:ilen(patname))
2394       if (lentmp.ne.0)
2395      &write (iout,*)" DIRTMP                          : ",
2396      &  tmpdir(:lentmp)
2397       write (iout,'(80(1h-))')
2398       endif
2399       return
2400       end
2401 c----------------------------------------------------------------------------
2402       subroutine card_concat(card)
2403       implicit real*8 (a-h,o-z)
2404       include 'DIMENSIONS'
2405       include 'COMMON.IOUNITS'
2406       character*(*) card
2407       character*80 karta,ucase
2408       external ilen
2409       read (inp,'(a)') karta
2410       karta=ucase(karta)
2411       card=' '
2412       do while (karta(80:80).eq.'&')
2413         card=card(:ilen(card)+1)//karta(:79)
2414         read (inp,'(a)') karta
2415         karta=ucase(karta)
2416       enddo
2417       card=card(:ilen(card)+1)//karta
2418       return
2419       end
2420 c----------------------------------------------------------------------------------
2421       subroutine readrst
2422       implicit real*8 (a-h,o-z)
2423       include 'DIMENSIONS'
2424       include 'COMMON.CHAIN'
2425       include 'COMMON.IOUNITS'
2426       include 'COMMON.MD'
2427       open(irest2,file=rest2name,status='unknown')
2428       read(irest2,*) totT,EK,potE,totE,t_bath
2429       totTafm=totT
2430       do i=1,2*nres
2431          read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2432       enddo
2433       do i=1,2*nres
2434          read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2435       enddo
2436       if(usampl) then
2437              read (irest2,*) iset
2438       endif
2439       close(irest2)
2440       return
2441       end
2442 c---------------------------------------------------------------------------------
2443       subroutine read_fragments
2444       implicit real*8 (a-h,o-z)
2445       include 'DIMENSIONS'
2446 #ifdef MPI
2447       include 'mpif.h'
2448 #endif
2449       include 'COMMON.SETUP'
2450       include 'COMMON.CHAIN'
2451       include 'COMMON.IOUNITS'
2452       include 'COMMON.MD'
2453       include 'COMMON.CONTROL'
2454       read(inp,*) nset,nfrag,npair,nfrag_back
2455       if(me.eq.king.or..not.out1file)
2456      & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2457      &  " nfrag_back",nfrag_back
2458       do iset=1,nset
2459          read(inp,*) mset(iset)
2460        do i=1,nfrag
2461          read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2462      &     qinfrag(i,iset)
2463          if(me.eq.king.or..not.out1file)
2464      &    write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2465      &     ifrag(2,i,iset), qinfrag(i,iset)
2466        enddo
2467        do i=1,npair
2468         read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2469      &    qinpair(i,iset)
2470         if(me.eq.king.or..not.out1file)
2471      &   write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2472      &    ipair(2,i,iset), qinpair(i,iset)
2473        enddo
2474        do i=1,nfrag_back
2475         read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2476      &     wfrag_back(3,i,iset),
2477      &     ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2478         if(me.eq.king.or..not.out1file)
2479      &   write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2480      &   wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2481        enddo
2482       enddo
2483       return
2484       end
2485 C---------------------------------------------------------------------------
2486       subroutine read_afminp
2487             implicit real*8 (a-h,o-z)
2488       include 'DIMENSIONS'
2489 #ifdef MPI
2490       include 'mpif.h'
2491 #endif
2492       include 'COMMON.SETUP'
2493       include 'COMMON.CONTROL'
2494       include 'COMMON.CHAIN'
2495       include 'COMMON.IOUNITS'
2496       include 'COMMON.SBRIDGE'
2497       character*320 afmcard
2498       print *, "wchodze"
2499       call card_concat(afmcard)
2500       call readi(afmcard,"BEG",afmbeg,0)
2501       call readi(afmcard,"END",afmend,0)
2502       call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
2503       call reada(afmcard,"VEL",velAFMconst,0.0d0)
2504       print *,'FORCE=' ,forceAFMconst
2505 CCCC NOW PROPERTIES FOR AFM
2506        distafminit=0.0d0
2507        do i=1,3
2508         distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
2509        enddo
2510         distafminit=dsqrt(distafminit)
2511         print *,'initdist',distafminit
2512       return
2513       end
2514
2515 c-------------------------------------------------------------------------------
2516       subroutine read_dist_constr
2517       implicit real*8 (a-h,o-z)
2518       include 'DIMENSIONS'
2519 #ifdef MPI
2520       include 'mpif.h'
2521 #endif
2522       include 'COMMON.SETUP'
2523       include 'COMMON.CONTROL'
2524       include 'COMMON.CHAIN'
2525       include 'COMMON.IOUNITS'
2526       include 'COMMON.SBRIDGE'
2527       integer ifrag_(2,100),ipair_(2,100)
2528       double precision wfrag_(100),wpair_(100)
2529       character*500 controlcard
2530       print *, "WCHODZE"
2531       write (iout,*) "Calling read_dist_constr"
2532 c      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2533 c      call flush(iout)
2534       if ((genconstr.gt.0).and.(constr_dist.eq.11)) then
2535       call gen_dist_constr
2536       go to 1712
2537       endif
2538       call card_concat(controlcard)
2539       call readi(controlcard,"NFRAG",nfrag_,0)
2540       call readi(controlcard,"NPAIR",npair_,0)
2541       call readi(controlcard,"NDIST",ndist_,0)
2542       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2543       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2544       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2545       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2546       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2547 c      write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2548 c      write (iout,*) "IFRAG"
2549 c      do i=1,nfrag_
2550 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2551 c      enddo
2552 c      write (iout,*) "IPAIR"
2553 c      do i=1,npair_
2554 c        write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2555 c      enddo
2556       call flush(iout)
2557       do i=1,nfrag_
2558         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2559         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2560      &    ifrag_(2,i)=nstart_sup+nsup-1
2561 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2562         call flush(iout)
2563         if (wfrag_(i).gt.0.0d0) then
2564         do j=ifrag_(1,i),ifrag_(2,i)-1
2565           do k=j+1,ifrag_(2,i)
2566 c            write (iout,*) "j",j," k",k
2567             ddjk=dist(j,k)
2568             if (constr_dist.eq.1) then
2569             nhpb=nhpb+1
2570             ihpb(nhpb)=j
2571             jhpb(nhpb)=k
2572               dhpb(nhpb)=ddjk
2573             forcon(nhpb)=wfrag_(i)
2574             else if (constr_dist.eq.2) then
2575               if (ddjk.le.dist_cut) then
2576                 nhpb=nhpb+1
2577                 ihpb(nhpb)=j
2578                 jhpb(nhpb)=k
2579                 dhpb(nhpb)=ddjk
2580                 forcon(nhpb)=wfrag_(i)
2581               endif
2582             else
2583               nhpb=nhpb+1
2584               ihpb(nhpb)=j
2585               jhpb(nhpb)=k
2586               dhpb(nhpb)=ddjk
2587               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2588             endif
2589 #ifdef MPI
2590             if (.not.out1file .or. me.eq.king)
2591      &      write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2592      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2593 #else
2594             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2595      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2596 #endif
2597           enddo
2598         enddo
2599         endif
2600       enddo
2601       do i=1,npair_
2602         if (wpair_(i).gt.0.0d0) then
2603         ii = ipair_(1,i)
2604         jj = ipair_(2,i)
2605         if (ii.gt.jj) then
2606           itemp=ii
2607           ii=jj
2608           jj=itemp
2609         endif
2610         do j=ifrag_(1,ii),ifrag_(2,ii)
2611           do k=ifrag_(1,jj),ifrag_(2,jj)
2612             nhpb=nhpb+1
2613             ihpb(nhpb)=j
2614             jhpb(nhpb)=k
2615             forcon(nhpb)=wpair_(i)
2616             dhpb(nhpb)=dist(j,k)
2617 #ifdef MPI
2618             if (.not.out1file .or. me.eq.king)
2619      &      write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2620      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2621 #else
2622             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2623      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2624 #endif
2625           enddo
2626         enddo
2627         endif
2628       enddo
2629       print *,ndist_
2630       do i=1,ndist_
2631         if (constr_dist.eq.11) then
2632         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2633      &     ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
2634         fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
2635         else
2636 C        print *,"in else"
2637         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2638      &     ibecarb(i),forcon(nhpb+1)
2639         endif
2640         if (forcon(nhpb+1).gt.0.0d0) then
2641           nhpb=nhpb+1
2642           if (ibecarb(i).gt.0) then
2643             ihpb(i)=ihpb(i)+nres
2644             jhpb(i)=jhpb(i)+nres
2645           endif
2646           if (dhpb(nhpb).eq.0.0d0)
2647      &       dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2648         endif
2649 C        read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2650 C        if (forcon(nhpb+1).gt.0.0d0) then
2651 C          nhpb=nhpb+1
2652 C          dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2653 #ifdef MPI
2654           if (.not.out1file .or. me.eq.king)
2655      &    write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2656      &     nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2657 #else
2658           write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2659      &     nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2660 #endif
2661
2662       enddo
2663  1712 continue
2664       call flush(iout)
2665       return
2666       end
2667 c-------------------------------------------------------------------------------
2668 #ifdef WINIFL
2669       subroutine flush(iu)
2670       return
2671       end
2672 #endif
2673 #ifdef AIX
2674       subroutine flush(iu)
2675       call flush_(iu)
2676       return
2677       end
2678 #endif
2679 c------------------------------------------------------------------------------
2680       subroutine copy_to_tmp(source)
2681       include "DIMENSIONS"
2682       include "COMMON.IOUNITS"
2683       character*(*) source
2684       character* 256 tmpfile
2685       integer ilen
2686       external ilen
2687       logical ex
2688       tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
2689       inquire(file=tmpfile,exist=ex)
2690       if (ex) then
2691         write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
2692      &   " to temporary directory..."
2693         write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
2694         call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
2695       endif
2696       return
2697       end
2698 c------------------------------------------------------------------------------
2699       subroutine move_from_tmp(source)
2700       include "DIMENSIONS"
2701       include "COMMON.IOUNITS"
2702       character*(*) source
2703       integer ilen
2704       external ilen
2705       write (*,*) "Moving ",source(:ilen(source)),
2706      & " from temporary directory to working directory"
2707       write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
2708       call system("/bin/mv "//source(:ilen(source))//" "//curdir)
2709       return
2710       end
2711 c------------------------------------------------------------------------------
2712       subroutine random_init(seed)
2713 C
2714 C Initialize random number generator
2715 C
2716       implicit real*8 (a-h,o-z)
2717       include 'DIMENSIONS'
2718 #ifdef MPI
2719       include 'mpif.h'
2720       logical OKRandom, prng_restart
2721       real*8  r1
2722       integer iseed_array(4)
2723 #endif
2724       include 'COMMON.IOUNITS'
2725       include 'COMMON.TIME1'
2726       include 'COMMON.THREAD'
2727       include 'COMMON.SBRIDGE'
2728       include 'COMMON.CONTROL'
2729       include 'COMMON.MCM'
2730       include 'COMMON.MAP'
2731       include 'COMMON.HEADER'
2732       include 'COMMON.CSA'
2733       include 'COMMON.CHAIN'
2734       include 'COMMON.MUCA'
2735       include 'COMMON.MD'
2736       include 'COMMON.FFIELD'
2737       include 'COMMON.SETUP'
2738       iseed=-dint(dabs(seed))
2739       if (iseed.eq.0) then
2740         write (iout,'(/80(1h*)/20x,a/80(1h*))')
2741      &    'Random seed undefined. The program will stop.'
2742         write (*,'(/80(1h*)/20x,a/80(1h*))')
2743      &    'Random seed undefined. The program will stop.'
2744 #ifdef MPI
2745         call mpi_finalize(mpi_comm_world,ierr)
2746 #endif
2747         stop 'Bad random seed.'
2748       endif
2749 #ifdef MPI
2750       if (fg_rank.eq.0) then
2751       seed=seed*(me+1)+1
2752 #ifdef AMD64
2753       if(me.eq.king)
2754      &  write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
2755       OKRandom = prng_restart(me,iseed)
2756 #else
2757       do i=1,4
2758        tmp=65536.0d0**(4-i)
2759        iseed_array(i) = dint(seed/tmp)
2760        seed=seed-iseed_array(i)*tmp
2761       enddo
2762       if(me.eq.king)
2763      & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
2764      &                 (iseed_array(i),i=1,4)
2765       write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
2766      &                 (iseed_array(i),i=1,4)
2767       OKRandom = prng_restart(me,iseed_array)
2768 #endif
2769       if (OKRandom) then
2770 c        r1 = prng_next(me)
2771         r1=ran_number(0.0D0,1.0D0)
2772         if(me.eq.king)
2773      &   write (iout,*) 'ran_num',r1
2774         if (r1.lt.0.0d0) OKRandom=.false.
2775       endif
2776       if (.not.OKRandom) then
2777         write (iout,*) 'PRNG IS NOT WORKING!!!'
2778         print *,'PRNG IS NOT WORKING!!!'
2779         if (me.eq.0) then
2780          call flush(iout)
2781          call mpi_abort(mpi_comm_world,error_msg,ierr)
2782          stop
2783         else
2784          write (iout,*) 'too many processors for parallel prng'
2785          write (*,*) 'too many processors for parallel prng'
2786          call flush(iout)
2787          stop
2788         endif
2789       endif
2790       endif
2791 #else
2792       call vrndst(iseed)
2793       write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
2794 #endif
2795       return
2796       end