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