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