source/unres/src_CSA_DiL/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       logical file_exist
  12 C Read force-field parameters except weights
  13       call parmread
  14 C Read job setup parameters
  15       call read_control
  16 C Read control parameters for energy minimzation if required
  17       if (minim) call read_minim
  18 C Read MCM control parameters if required
  19 c      if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
  20 C Read MD control parameters if reqjuired
  21 c      if (modecalc.eq.12) call read_MDpar
  22 C Read MREMD control parameters if required
  23 c      if (modecalc.eq.14) then
  24 c         call read_MDpar
  25 c         call read_REMDpar
  26 c      endif
  27 C Read MUCA control parameters if required
  28 c      if (lmuca) call read_muca
  29 C Read CSA control parameters if required (from fort.40 if exists
  30 C otherwise from general input file)
  31       if (modecalc.eq.8) then
  32        inquire (file="fort.40",exist=file_exist)
  33        if (.not.file_exist) call csaread
  34       endif
  35 cfmc      if (modecalc.eq.10) call mcmfread
  36 C Read molecule information, molecule geometry, energy-term weights, and
  37 C restraints if requested
  38       call molread
  39
  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 c      include 'COMMON.THREAD'
  71       include 'COMMON.SBRIDGE'
  72       include 'COMMON.CONTROL'
  73       include 'COMMON.MCM'
  74 c      include 'COMMON.MAP'
  75       include 'COMMON.HEADER'
  76 c      include 'COMMON.CSA'
  77       include 'COMMON.CHAIN'
  78 c      include 'COMMON.MUCA'
  79 c      include 'COMMON.MD'
  80       include 'COMMON.FFIELD'
  81       include 'COMMON.SETUP'
  82       COMMON /MACHSW/ KDIAG,ICORFL,IXDR
  83       character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
  84       character*80 ucase
  85       character*320 controlcard
  86
  87       nglob_csa=0
  88       eglob_csa=1d99
  89       nmin_csa=0
  90       read (INP,'(a)') titel
  91       call card_concat(controlcard)
  92 c      out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
  93 c      print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
  94       call reada(controlcard,'SEED',seed,0.0D0)
  95       call random_init(seed)
  96 C Set up the time limit (caution! The time must be input in minutes!)
  97       read_cart=index(controlcard,'READ_CART').gt.0
  98       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
  99       call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
 100       unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
 101       call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
 102       call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
 103       call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
 104       call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
 105       call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
 106       call reada(controlcard,'DRMS',drms,0.1D0)
 107       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 108        write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
 109        write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
 110        write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
 111        write (iout,'(a,f10.1)')'DRMS    = ',drms
 112        write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
 113        write (iout,'(a,f10.1)') 'Time limit (min):',timlim
 114       endif
 115       call readi(controlcard,'NZ_START',nz_start,0)
 116       call readi(controlcard,'NZ_END',nz_end,0)
 117       call readi(controlcard,'IZ_SC',iz_sc,0)
 118       timlim=60.0D0*timlim
 119       safety = 60.0d0*safety
 120       timem=timlim
 121       modecalc=0
 122       call reada(controlcard,"T_BATH",t_bath,300.0d0)
 123       minim=(index(controlcard,'MINIMIZE').gt.0)
 124       dccart=(index(controlcard,'CART').gt.0)
 125       overlapsc=(index(controlcard,'OVERLAP').gt.0)
 126       overlapsc=.not.overlapsc
 127       searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
 128       searchsc=.not.searchsc
 129       sideadd=(index(controlcard,'SIDEADD').gt.0)
 130       energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
 131       outpdb=(index(controlcard,'PDBOUT').gt.0)
 132       outmol2=(index(controlcard,'MOL2OUT').gt.0)
 133       pdbref=(index(controlcard,'PDBREF').gt.0)
 134       refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
 135       indpdb=index(controlcard,'PDBSTART')
 136       extconf=(index(controlcard,'EXTCONF').gt.0)
 137       call readi(controlcard,'IPRINT',iprint,0)
 138       call readi(controlcard,'MAXGEN',maxgen,10000)
 139       call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
 140       call readi(controlcard,"KDIAG",kdiag,0)
 141       call readi(controlcard,"RESCALE_MODE",rescale_mode,0)
 142       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
 143      & write (iout,*) "RESCALE_MODE",rescale_mode
 144       split_ene=index(controlcard,'SPLIT_ENE').gt.0
 145       if (index(controlcard,'REGULAR').gt.0.0D0) then
 146         call reada(controlcard,'WEIDIS',weidis,0.1D0)
 147         modecalc=1
 148         refstr=.true.
 149       endif
 150       if (index(controlcard,'CHECKGRAD').gt.0) then
 151         modecalc=5
 152         if (index(controlcard,'CART').gt.0) then
 153           icheckgrad=1
 154         elseif (index(controlcard,'CARINT').gt.0) then
 155           icheckgrad=2
 156         else
 157           icheckgrad=3
 158         endif
 159       elseif (index(controlcard,'THREAD').gt.0) then
 160         modecalc=2
 161         call readi(controlcard,'THREAD',nthread,0)
 162         if (nthread.gt.0) then
 163           call reada(controlcard,'WEIDIS',weidis,0.1D0)
 164         else
 165           if (fg_rank.eq.0)
 166      &    write (iout,'(a)')'A number has to follow the THREAD keyword.'
 167           stop 'Error termination in Read_Control.'
 168         endif
 169       else if (index(controlcard,'MCMA').gt.0) then
 170         modecalc=3
 171       else if (index(controlcard,'MCEE').gt.0) then
 172         modecalc=6
 173       else if (index(controlcard,'MULTCONF').gt.0) then
 174         modecalc=4
 175       else if (index(controlcard,'MAP').gt.0) then
 176         modecalc=7
 177         call readi(controlcard,'MAP',nmap,0)
 178       else if (index(controlcard,'CSA').gt.0) then
 179         modecalc=8
 180 crc      else if (index(controlcard,'ZSCORE').gt.0) then
 181 crc
 182 crc  ZSCORE is rm from UNRES, modecalc=9 is available
 183 crc
 184 crc        modecalc=9
 185 cfcm      else if (index(controlcard,'MCMF').gt.0) then
 186 cfmc        modecalc=10
 187       else if (index(controlcard,'SOFTREG').gt.0) then
 188         modecalc=11
 189       else if (index(controlcard,'CHECK_BOND').gt.0) then
 190         modecalc=-1
 191       else if (index(controlcard,'TEST').gt.0) then
 192         modecalc=-2
 193       else if (index(controlcard,'MD').gt.0) then
 194         modecalc=12
 195       else if (index(controlcard,'RE ').gt.0) then
 196         modecalc=14
 197       endif
 198
 199       lmuca=index(controlcard,'MUCA').gt.0
 200       call readi(controlcard,'MUCADYN',mucadyn,0)
 201       call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
 202       if (lmuca .and. (me.eq.king .or. .not.out1file ))
 203      & then
 204        write (iout,*) 'MUCADYN=',mucadyn
 205        write (iout,*) 'MUCASMOOTH=',muca_smooth
 206       endif
 207
 208       iscode=index(controlcard,'ONE_LETTER')
 209       indphi=index(controlcard,'PHI')
 210       indback=index(controlcard,'BACK')
 211       iranconf=index(controlcard,'RAND_CONF')
 212       i2ndstr=index(controlcard,'USE_SEC_PRED')
 213       gradout=index(controlcard,'GRADOUT').gt.0
 214       gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
 215
 216       if(me.eq.king.or..not.out1file)
 217      & write (iout,'(2a)') diagmeth(kdiag),
 218      &  ' routine used to diagonalize matrices.'
 219       return
 220       end
 221 c------------------------------------------------------------------------------
 222       subroutine molread
 223 C
 224 C Read molecular data.
 225 C
 226       implicit real*8 (a-h,o-z)
 227       include 'DIMENSIONS'
 228 #ifdef MPI
 229       include 'mpif.h'
 230       integer error_msg
 231 #endif
 232       include 'COMMON.IOUNITS'
 233       include 'COMMON.GEO'
 234       include 'COMMON.VAR'
 235       include 'COMMON.INTERACT'
 236       include 'COMMON.LOCAL'
 237       include 'COMMON.NAMES'
 238       include 'COMMON.CHAIN'
 239       include 'COMMON.FFIELD'
 240       include 'COMMON.SBRIDGE'
 241       include 'COMMON.HEADER'
 242       include 'COMMON.CONTROL'
 243 c      include 'COMMON.DBASE'
 244 c      include 'COMMON.THREAD'
 245       include 'COMMON.CONTACTS'
 246       include 'COMMON.TORCNSTR'
 247       include 'COMMON.TIME1'
 248       include 'COMMON.BOUNDS'
 249 c      include 'COMMON.MD'
 250 c      include 'COMMON.REMD'
 251       include 'COMMON.SETUP'
 252       character*4 sequence(maxres)
 253       integer rescode
 254       double precision x(maxvar)
 255       character*256 pdbfile
 256       character*320 weightcard
 257       character*80 weightcard_t,ucase
 258       dimension itype_pdb(maxres)
 259       common /pizda/ itype_pdb
 260       logical seq_comp,fail
 261       double precision energia(0:n_ene)
 262       integer ilen
 263       external ilen
 264 C
 265 C Body
 266 C
 267 C Read weights of the subsequent energy terms.
 268
 269        call card_concat(weightcard)
 270        call reada(weightcard,'WLONG',wlong,1.0D0)
 271        call reada(weightcard,'WSC',wsc,wlong)
 272        call reada(weightcard,'WSCP',wscp,wlong)
 273        call reada(weightcard,'WELEC',welec,1.0D0)
 274        call reada(weightcard,'WVDWPP',wvdwpp,welec)
 275        call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
 276        call reada(weightcard,'WCORR4',wcorr4,0.0D0)
 277        call reada(weightcard,'WCORR5',wcorr5,0.0D0)
 278        call reada(weightcard,'WCORR6',wcorr6,0.0D0)
 279        call reada(weightcard,'WTURN3',wturn3,1.0D0)
 280        call reada(weightcard,'WTURN4',wturn4,1.0D0)
 281        call reada(weightcard,'WTURN6',wturn6,1.0D0)
 282        call reada(weightcard,'WSCCOR',wsccor,1.0D0)
 283        call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
 284        call reada(weightcard,'WBOND',wbond,1.0D0)
 285        call reada(weightcard,'WTOR',wtor,1.0D0)
 286        call reada(weightcard,'WTORD',wtor_d,1.0D0)
 287        call reada(weightcard,'WANG',wang,1.0D0)
 288        call reada(weightcard,'WSCLOC',wscloc,1.0D0)
 289 C     Juyong
 290        call reada(weightcard,'WDFAD',wdfa_dist,0.0d0)
 291        call reada(weightcard,'WDFAT',wdfa_tor,0.0d0)
 292        call reada(weightcard,'WDFAN',wdfa_nei,0.0d0)
 293        call reada(weightcard,'WDFAB',wdfa_beta,0.0d0)
 294 C
 295        call reada(weightcard,'SCAL14',scal14,0.4D0)
 296        call reada(weightcard,'SCALSCP',scalscp,1.0d0)
 297        call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
 298        call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
 299        call reada(weightcard,'TEMP0',temp0,300.0d0)
 300        if (index(weightcard,'SOFT').gt.0) ipot=6
 301 C 12/1/95 Added weight for the multi-body term WCORR
 302        call reada(weightcard,'WCORRH',wcorr,1.0D0)
 303        if (wcorr4.gt.0.0d0) wcorr=wcorr4
 304        weights(1)=wsc
 305        weights(2)=wscp
 306        weights(3)=welec
 307        weights(4)=wcorr
 308        weights(5)=wcorr5
 309        weights(6)=wcorr6
 310        weights(7)=wel_loc
 311        weights(8)=wturn3
 312        weights(9)=wturn4
 313        weights(10)=wturn6
 314        weights(11)=wang
 315        weights(12)=wscloc
 316        weights(13)=wtor
 317        weights(14)=wtor_d
 318        weights(15)=wstrain
 319        weights(16)=wvdwpp
 320        weights(17)=wbond
 321        weights(18)=scal14
 322        weights(21)=wsccor
 323 C     JUYONG
 324        weights(24)=wdfa_dist
 325        weights(25)=wdfa_tor
 326        weights(26)=wdfa_nei
 327        weights(27)=wdfa_beta
 328 C
 329
 330       if(me.eq.king.or..not.out1file)
 331      & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
 332      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
 333      &  wturn4,wturn6,
 334      &  wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta
 335
 336    10 format (/'Energy-term weights (unscaled):'//
 337      & 'WSCC=   ',f10.6,' (SC-SC)'/
 338      & 'WSCP=   ',f10.6,' (SC-p)'/
 339      & 'WELEC=  ',f10.6,' (p-p electr)'/
 340      & 'WVDWPP= ',f10.6,' (p-p VDW)'/
 341      & 'WBOND=  ',f10.6,' (stretching)'/
 342      & 'WANG=   ',f10.6,' (bending)'/
 343      & 'WSCLOC= ',f10.6,' (SC local)'/
 344      & 'WTOR=   ',f10.6,' (torsional)'/
 345      & 'WTORD=  ',f10.6,' (double torsional)'/
 346      & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
 347      & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
 348      & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
 349      & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
 350      & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
 351      & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
 352      & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
 353      & 'WTURN4= ',f10.6,' (turns, 4th order)'/
 354      & 'WTURN6= ',f10.6,' (turns, 6th order)'/
 355      & 'WDFA_D= ',f10.6,' (DFA, distance)'   /
 356      & 'WDFA_T= ',f10.6,' (DFA, torsional)'   /
 357      & 'WDFA_N= ',f10.6,' (DFA, number of neighbor)'   /
 358      & 'WDFA_B= ',f10.6,' (DFA, beta formation)')
 359
 360       if(me.eq.king.or..not.out1file)then
 361        if (wcorr4.gt.0.0d0) then
 362         write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
 363      &   'between contact pairs of peptide groups'
 364         write (iout,'(2(a,f5.3/))')
 365      &  'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
 366      &  'Range of quenching the correlation terms:',2*delt_corr
 367        else if (wcorr.gt.0.0d0) then
 368         write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
 369      &   'between contact pairs of peptide groups'
 370        endif
 371        write (iout,'(a,f8.3)')
 372      &  'Scaling factor of 1,4 SC-p interactions:',scal14
 373        write (iout,'(a,f8.3)')
 374      &  'General scaling factor of SC-p interactions:',scalscp
 375       endif
 376       r0_corr=cutoff_corr-delt_corr
 377       do i=1,20
 378         aad(i,1)=scalscp*aad(i,1)
 379         aad(i,2)=scalscp*aad(i,2)
 380         bad(i,1)=scalscp*bad(i,1)
 381         bad(i,2)=scalscp*bad(i,2)
 382       enddo
 383 c      call rescale_weights(t_bath)
 384       if(me.eq.king.or..not.out1file)
 385      & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
 386      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
 387      &  wturn4,wturn6,
 388      &  wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta
 389
 390    22 format (/'Energy-term weights (scaled):'//
 391      & 'WSCC=   ',f10.6,' (SC-SC)'/
 392      & 'WSCP=   ',f10.6,' (SC-p)'/
 393      & 'WELEC=  ',f10.6,' (p-p electr)'/
 394      & 'WVDWPP= ',f10.6,' (p-p VDW)'/
 395      & 'WBOND=  ',f10.6,' (stretching)'/
 396      & 'WANG=   ',f10.6,' (bending)'/
 397      & 'WSCLOC= ',f10.6,' (SC local)'/
 398      & 'WTOR=   ',f10.6,' (torsional)'/
 399      & 'WTORD=  ',f10.6,' (double torsional)'/
 400      & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
 401      & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
 402      & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
 403      & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
 404      & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
 405      & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
 406      & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
 407      & 'WTURN4= ',f10.6,' (turns, 4th order)'/
 408      & 'WTURN6= ',f10.6,' (turns, 6th order)'/
 409      & 'WDFA_D= ',f10.6,' (DFA, distance)'   /
 410      & 'WDFA_T= ',f10.6,' (DFA, torsional)'   /
 411      & 'WDFA_N= ',f10.6,' (DFA, number of neighbor)'   /
 412      & 'WDFA_B= ',f10.6,' (DFA, beta formation)')
 413
 414       if(me.eq.king.or..not.out1file)
 415      & write (iout,*) "Reference temperature for weights calculation:",
 416      &  temp0
 417       call reada(weightcard,"D0CM",d0cm,3.78d0)
 418       call reada(weightcard,"AKCM",akcm,15.1d0)
 419       call reada(weightcard,"AKTH",akth,11.0d0)
 420       call reada(weightcard,"AKCT",akct,12.0d0)
 421       call reada(weightcard,"V1SS",v1ss,-1.08d0)
 422       call reada(weightcard,"V2SS",v2ss,7.61d0)
 423       call reada(weightcard,"V3SS",v3ss,13.7d0)
 424       call reada(weightcard,"EBR",ebr,-5.50D0)
 425       if(me.eq.king.or..not.out1file) then
 426        write (iout,*) "Parameters of the SS-bond potential:"
 427        write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
 428      & " AKCT",akct
 429        write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
 430        write (iout,*) "EBR",ebr
 431        print *,'indpdb=',indpdb,' pdbref=',pdbref
 432       endif
 433       if (indpdb.gt.0 .or. pdbref) then
 434         read(inp,'(a)') pdbfile
 435         if(me.eq.king.or..not.out1file)
 436      &   write (iout,'(2a)') 'PDB data will be read from file ',
 437      &   pdbfile(:ilen(pdbfile))
 438         open(ipdbin,file=pdbfile,status='old',err=33)
 439         goto 34
 440   33    write (iout,'(a)') 'Error opening PDB file.'
 441         stop
 442   34    continue
 443 c        print *,'Begin reading pdb data'
 444         call readpdb
 445 c        print *,'Finished reading pdb data'
 446         if(me.eq.king.or..not.out1file)
 447      &   write (iout,'(a,i3,a,i3)')'nsup=',nsup,
 448      &   ' nstart_sup=',nstart_sup
 449         do i=1,nres
 450           itype_pdb(i)=itype(i)
 451         enddo
 452         close (ipdbin)
 453         nnt=nstart_sup
 454         nct=nstart_sup+nsup-1
 455         call contact(.false.,ncont_ref,icont_ref,co)
 456
 457         if (sideadd) then
 458          if(me.eq.king.or..not.out1file)
 459      &    write(iout,*)'Adding sidechains'
 460          maxsi=1000
 461          do i=2,nres-1
 462           iti=itype(i)
 463           if (iti.ne.10) then
 464             nsi=0
 465             fail=.true.
 466             do while (fail.and.nsi.le.maxsi)
 467 c              call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
 468               nsi=nsi+1
 469             enddo
 470             if(fail) write(iout,*)'Adding sidechain failed for res ',
 471      &              i,' after ',nsi,' trials'
 472           endif
 473          enddo
 474         endif
 475       endif
 476
 477       if (indpdb.eq.0) then
 478 C Read sequence if not taken from the pdb file.
 479         read (inp,*) nres
 480 c        print *,'nres=',nres
 481         if (iscode.gt.0) then
 482           read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
 483         else
 484           read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
 485         endif
 486 C Convert sequence to numeric code
 487         do i=1,nres
 488           itype(i)=rescode(i,sequence(i),iscode)
 489         enddo
 490 C Assign initial virtual bond lengths
 491         do i=2,nres
 492           vbld(i)=vbl
 493           vbld_inv(i)=vblinv
 494         enddo
 495         do i=2,nres-1
 496           vbld(i+nres)=dsc(iabs(itype(i)))
 497           vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
 498 c          write (iout,*) "i",i," itype",itype(i),
 499 c     &      " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
 500         enddo
 501       endif
 502 c      print *,nres
 503 c      print '(20i4)',(itype(i),i=1,nres)
 504       do i=1,nres
 505 #ifdef PROCOR
 506         if (itype(i).eq.21 .or. itype(i+1).eq.21) then
 507 #else
 508         if (itype(i).eq.21) then
 509 #endif
 510           itel(i)=0
 511 #ifdef PROCOR
 512         else if (iabs(itype(i+1)).ne.20) then
 513 #else
 514         else if (iabs(itype(i)).ne.20) then
 515 #endif
 516           itel(i)=1
 517         else
 518           itel(i)=2
 519         endif
 520       enddo
 521       if(me.eq.king.or..not.out1file)then
 522        write (iout,*) "ITEL"
 523        do i=1,nres-1
 524          write (iout,*) i,itype(i),itel(i)
 525        enddo
 526        print *,'Call Read_Bridge.'
 527       endif
 528       call read_bridge
 529 C 8/13/98 Set limits to generating the dihedral angles
 530       do i=1,nres
 531         phibound(1,i)=-pi
 532         phibound(2,i)=pi
 533       enddo
 534       read (inp,*) ndih_constr
 535       if (ndih_constr.gt.0) then
 536         read (inp,*) ftors
 537         read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
 538         if(me.eq.king.or..not.out1file)then
 539          write (iout,*)
 540      &   'There are',ndih_constr,' constraints on phi angles.'
 541          do i=1,ndih_constr
 542           write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
 543          enddo
 544         endif
 545         do i=1,ndih_constr
 546           phi0(i)=deg2rad*phi0(i)
 547           drange(i)=deg2rad*drange(i)
 548         enddo
 549         if(me.eq.king.or..not.out1file)
 550      &   write (iout,*) 'FTORS',ftors
 551         do i=1,ndih_constr
 552           ii = idih_constr(i)
 553           phibound(1,ii) = phi0(i)-drange(i)
 554           phibound(2,ii) = phi0(i)+drange(i)
 555         enddo
 556       endif
 557       nnt=1
 558 #ifdef MPI
 559       if (me.eq.king) then
 560 #endif
 561        write (iout,'(a)') 'Boundaries in phi angle sampling:'
 562        do i=1,nres
 563          write (iout,'(a3,i5,2f10.1)')
 564      &   restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
 565        enddo
 566 #ifdef MP
 567       endif
 568 #endif
 569       nct=nres
 570 cd      print *,'NNT=',NNT,' NCT=',NCT
 571       if (itype(1).eq.21) nnt=2
 572       if (itype(nres).eq.21) nct=nct-1
 573
 574 C     Juyong:READ init_vars
 575 C     Initialize variables!
 576 C     Juyong:READ read_info
 577 C     READ fragment information!!
 578 C     both routines should be in dfa.F file!!
 579
 580       if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
 581      &            wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
 582        call init_dfa_vars
 583        print*, 'init_dfa_vars finished!'
 584        call read_dfa_info
 585        print*, 'read_dfa_info finished!'
 586       endif
 587 C
 588 C
 589
 590
 591       if (pdbref) then
 592         if(me.eq.king.or..not.out1file)
 593      &   write (iout,'(a,i3)') 'nsup=',nsup
 594         nstart_seq=nnt
 595         if (nsup.le.(nct-nnt+1)) then
 596           do i=0,nct-nnt+1-nsup
 597             if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
 598               nstart_seq=nnt+i
 599               goto 111
 600             endif
 601           enddo
 602           write (iout,'(a)')
 603      &            'Error - sequences to be superposed do not match.'
 604           stop
 605         else
 606           do i=0,nsup-(nct-nnt+1)
 607             if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
 608      &      then
 609               nstart_sup=nstart_sup+i
 610               nsup=nct-nnt+1
 611               goto 111
 612             endif
 613           enddo
 614           write (iout,'(a)')
 615      &            'Error - sequences to be superposed do not match.'
 616         endif
 617   111   continue
 618         if (nsup.eq.0) nsup=nct-nnt
 619         if (nstart_sup.eq.0) nstart_sup=nnt
 620         if (nstart_seq.eq.0) nstart_seq=nnt
 621         if(me.eq.king.or..not.out1file)
 622      &   write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
 623      &                 ' nstart_seq=',nstart_seq
 624       endif
 625 c--- Zscore rms -------
 626       if (nz_start.eq.0) nz_start=nnt
 627       if (nz_end.eq.0 .and. nsup.gt.0) then
 628         nz_end=nnt+nsup-1
 629       else if (nz_end.eq.0) then
 630         nz_end=nct
 631       endif
 632       if(me.eq.king.or..not.out1file)then
 633        write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
 634        write (iout,*) 'IZ_SC=',iz_sc
 635       endif
 636 c----------------------
 637       call init_int_table
 638       if (refstr) then
 639         if (.not.pdbref) then
 640           call read_angles(inp,*38)
 641           goto 39
 642    38     write (iout,'(a)') 'Error reading reference structure.'
 643 #ifdef MPI
 644           call MPI_Finalize(MPI_COMM_WORLD,IERROR)
 645           stop 'Error reading reference structure'
 646 #endif
 647    39     call chainbuild
 648           call setup_var
 649 czscore          call geom_to_var(nvar,coord_exp_zs(1,1))
 650           nstart_sup=nnt
 651           nstart_seq=nnt
 652           nsup=nct-nnt+1
 653           do i=1,2*nres
 654             do j=1,3
 655               cref(j,i)=c(j,i)
 656             enddo
 657           enddo
 658           call contact(.true.,ncont_ref,icont_ref,co)
 659         endif
 660 c        write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
 661         call flush(iout)
 662         if (constr_dist.gt.0) call read_dist_constr
 663 c        write (iout,*) "After read_dist_constr nhpb",nhpb
 664         call hpb_partition
 665         if(me.eq.king.or..not.out1file)
 666      &   write (iout,*) 'Contact order:',co
 667         if (pdbref) then
 668         if(me.eq.king.or..not.out1file)
 669      &   write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
 670         do i=1,ncont_ref
 671           do j=1,2
 672             icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
 673           enddo
 674           if(me.eq.king.or..not.out1file)
 675      &     write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
 676      &     icont_ref(1,i),' ',
 677      &     restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
 678         enddo
 679         endif
 680       endif
 681       if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
 682      &    .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
 683      &    modecalc.ne.10) then
 684 C If input structure hasn't been supplied from the PDB file read or generate
 685 C initial geometry.
 686         if (iranconf.eq.0 .and. .not. extconf) then
 687           if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
 688      &     write (iout,'(a)') 'Initial geometry will be read in.'
 689           if (read_cart) then
 690             read(inp,'(8f10.5)',end=36,err=36)
 691      &       ((c(l,k),l=1,3),k=1,nres),
 692      &       ((c(l,k+nres),l=1,3),k=nnt,nct)
 693             call int_from_cart1(.false.)
 694             do i=1,nres-1
 695               do j=1,3
 696                 dc(j,i)=c(j,i+1)-c(j,i)
 697                 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
 698               enddo
 699             enddo
 700             do i=nnt,nct
 701               if (itype(i).ne.10) then
 702                 do j=1,3
 703                   dc(j,i+nres)=c(j,i+nres)-c(j,i)
 704                   dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
 705                 enddo
 706               endif
 707             enddo
 708             return
 709           else
 710             call read_angles(inp,*36)
 711           endif
 712           goto 37
 713    36     write (iout,'(a)') 'Error reading angle file.'
 714 #ifdef MPI
 715           call mpi_finalize( MPI_COMM_WORLD,IERR )
 716 #endif
 717           stop 'Error reading angle file.'
 718    37     continue
 719         else if (extconf) then
 720          if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
 721      &    write (iout,'(a)') 'Extended chain initial geometry.'
 722          do i=3,nres
 723           theta(i)=90d0*deg2rad
 724          enddo
 725          do i=4,nres
 726           phi(i)=180d0*deg2rad
 727          enddo
 728          do i=2,nres-1
 729           alph(i)=110d0*deg2rad
 730          enddo
 731          do i=2,nres-1
 732           omeg(i)=-120d0*deg2rad
 733           if (itype(i).le.0) omeg(i)=-omeg(i)
 734          enddo
 735         else
 736           if(me.eq.king.or..not.out1file)
 737      &     write (iout,'(a)') 'Random-generated initial geometry.'
 738
 739
 740 #ifdef MPI
 741           if (me.eq.king  .or. fg_rank.eq.0 .and. (
 742      &           modecalc.eq.12 .or. modecalc.eq.14) ) then
 743 #endif
 744             do itrial=1,100
 745               itmp=1
 746 c              call gen_rand_conf(itmp,*30)
 747               goto 40
 748    30         write (iout,*) 'Failed to generate random conformation',
 749      &          ', itrial=',itrial
 750               write (*,*) 'Processor:',me,
 751      &          ' Failed to generate random conformation',
 752      &          ' itrial=',itrial
 753               call intout
 754
 755 #ifdef AIX
 756               call flush_(iout)
 757 #else
 758               call flush(iout)
 759 #endif
 760             enddo
 761             write (iout,'(a,i3,a)') 'Processor:',me,
 762      &        ' error in generating random conformation.'
 763             write (*,'(a,i3,a)') 'Processor:',me,
 764      &        ' error in generating random conformation.'
 765             call flush(iout)
 766 #ifdef MPI
 767             call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
 768    40       continue
 769           endif
 770 #else
 771           do itrial=1,100
 772             itmp=1
 773 c            call gen_rand_conf(itmp,*31)
 774             goto 40
 775    31       write (iout,*) 'Failed to generate random conformation',
 776      &        ', itrial=',itrial
 777             write (*,*) 'Failed to generate random conformation',
 778      &        ', itrial=',itrial
 779           enddo
 780           write (iout,'(a,i3,a)') 'Processor:',me,
 781      &      ' error in generating random conformation.'
 782           write (*,'(a,i3,a)') 'Processor:',me,
 783      &      ' error in generating random conformation.'
 784           stop
 785    40     continue
 786 #endif
 787         endif
 788       elseif (modecalc.eq.4) then
 789         read (inp,'(a)') intinname
 790         open (intin,file=intinname,status='old',err=333)
 791         if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
 792      &  write (iout,'(a)') 'intinname',intinname
 793         write (*,'(a)') 'Processor',myrank,' intinname',intinname
 794         goto 334
 795   333   write (iout,'(2a)') 'Error opening angle file ',intinname
 796 #ifdef MPI
 797         call MPI_Finalize(MPI_COMM_WORLD,IERR)
 798 #endif
 799         stop 'Error opening angle file.'
 800   334   continue
 801
 802       endif
 803 C Generate distance constraints, if the PDB structure is to be regularized.
 804       if (nthread.gt.0) then
 805         call read_threadbase
 806       endif
 807       call setup_var
 808       if (me.eq.king .or. .not. out1file)
 809      & call intout
 810       if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
 811         write (iout,'(/a,i3,a)')
 812      &  'The chain contains',ns,' disulfide-bridging cysteines.'
 813         write (iout,'(20i4)') (iss(i),i=1,ns)
 814         write (iout,'(/a/)') 'Pre-formed links are:'
 815         do i=1,nss
 816           i1=ihpb(i)-nres
 817           i2=jhpb(i)-nres
 818           it1=itype(i1)
 819           it2=itype(i2)
 820           if (me.eq.king.or..not.out1file)
 821      &    write (iout,'(2a,i3,3a,i3,a,3f10.3)')
 822      &    restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
 823      &    ebr,forcon(i)
 824         enddo
 825         write (iout,'(a)')
 826       endif
 827 c       if (i2ndstr.gt.0) call secstrp2dihc
 828 c      call geom_to_var(nvar,x)
 829 c      call etotal(energia(0))
 830 c      call enerprint(energia(0))
 831 c      call briefout(0,etot)
 832 c      stop
 833 cd    write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
 834 cd    write (iout,'(a)') 'Variable list:'
 835 cd    write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
 836 #ifdef MPI
 837       if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
 838      &  write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
 839      &  'Processor',myrank,': end reading molecular data.'
 840 #endif
 841       return
 842       end
 843 c--------------------------------------------------------------------------
 844       logical function seq_comp(itypea,itypeb,length)
 845       implicit none
 846       integer length,itypea(length),itypeb(length)
 847       integer i
 848       do i=1,length
 849         if (itypea(i).ne.itypeb(i)) then
 850           seq_comp=.false.
 851           return
 852         endif
 853       enddo
 854       seq_comp=.true.
 855       return
 856       end
 857 c-----------------------------------------------------------------------------
 858       subroutine read_bridge
 859 C Read information about disulfide bridges.
 860       implicit real*8 (a-h,o-z)
 861       include 'DIMENSIONS'
 862 #ifdef MPI
 863       include 'mpif.h'
 864 #endif
 865       include 'COMMON.IOUNITS'
 866       include 'COMMON.GEO'
 867       include 'COMMON.VAR'
 868       include 'COMMON.INTERACT'
 869       include 'COMMON.LOCAL'
 870       include 'COMMON.NAMES'
 871       include 'COMMON.CHAIN'
 872       include 'COMMON.FFIELD'
 873       include 'COMMON.SBRIDGE'
 874       include 'COMMON.HEADER'
 875       include 'COMMON.CONTROL'
 876 c      include 'COMMON.DBASE'
 877 c      include 'COMMON.THREAD'
 878       include 'COMMON.TIME1'
 879       include 'COMMON.SETUP'
 880 C Read bridging residues.
 881       read (inp,*) ns,(iss(i),i=1,ns)
 882       print *,'ns=',ns
 883       if(me.eq.king.or..not.out1file)
 884      &  write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
 885 C Check whether the specified bridging residues are cystines.
 886       do i=1,ns
 887         if (itype(iss(i)).ne.1) then
 888           if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
 889      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
 890      &   ' can form a disulfide bridge?!!!'
 891           write (*,'(2a,i3,a)')
 892      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
 893      &   ' can form a disulfide bridge?!!!'
 894 #ifdef MPI
 895          call MPI_Finalize(MPI_COMM_WORLD,ierror)
 896          stop
 897 #endif
 898         endif
 899       enddo
 900 C Read preformed bridges.
 901       if (ns.gt.0) then
 902       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
 903       write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
 904       if (nss.gt.0) then
 905         nhpb=nss
 906 C Check if the residues involved in bridges are in the specified list of
 907 C bridging residues.
 908         do i=1,nss
 909           do j=1,i-1
 910             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
 911      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
 912               write (iout,'(a,i3,a)') 'Disulfide pair',i,
 913      &      ' contains residues present in other pairs.'
 914               write (*,'(a,i3,a)') 'Disulfide pair',i,
 915      &      ' contains residues present in other pairs.'
 916 #ifdef MPI
 917               call MPI_Finalize(MPI_COMM_WORLD,ierror)
 918               stop
 919 #endif
 920             endif
 921           enddo
 922           do j=1,ns
 923             if (ihpb(i).eq.iss(j)) goto 10
 924           enddo
 925           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 926    10     continue
 927           do j=1,ns
 928             if (jhpb(i).eq.iss(j)) goto 20
 929           enddo
 930           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 931    20     continue
 932           dhpb(i)=dbr
 933           forcon(i)=fbr
 934         enddo
 935         do i=1,nss
 936           ihpb(i)=ihpb(i)+nres
 937           jhpb(i)=jhpb(i)+nres
 938         enddo
 939       endif
 940       endif
 941       return
 942       end
 943 c----------------------------------------------------------------------------
 944       subroutine read_x(kanal,*)
 945       implicit real*8 (a-h,o-z)
 946       include 'DIMENSIONS'
 947       include 'COMMON.GEO'
 948       include 'COMMON.VAR'
 949       include 'COMMON.CHAIN'
 950       include 'COMMON.IOUNITS'
 951       include 'COMMON.CONTROL'
 952       include 'COMMON.LOCAL'
 953       include 'COMMON.INTERACT'
 954 c Read coordinates from input
 955 c
 956       read(kanal,'(8f10.5)',end=10,err=10)
 957      &  ((c(l,k),l=1,3),k=1,nres),
 958      &  ((c(l,k+nres),l=1,3),k=nnt,nct)
 959       do j=1,3
 960         c(j,nres+1)=c(j,1)
 961         c(j,2*nres)=c(j,nres)
 962       enddo
 963       call int_from_cart1(.false.)
 964       do i=1,nres-1
 965         do j=1,3
 966           dc(j,i)=c(j,i+1)-c(j,i)
 967           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
 968         enddo
 969       enddo
 970       do i=nnt,nct
 971         if (itype(i).ne.10) then
 972           do j=1,3
 973             dc(j,i+nres)=c(j,i+nres)-c(j,i)
 974             dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
 975           enddo
 976         endif
 977       enddo
 978
 979       return
 980    10 return1
 981       end
 982 c------------------------------------------------------------------------------
 983       subroutine setup_var
 984       implicit real*8 (a-h,o-z)
 985       include 'DIMENSIONS'
 986       include 'COMMON.IOUNITS'
 987       include 'COMMON.GEO'
 988       include 'COMMON.VAR'
 989       include 'COMMON.INTERACT'
 990       include 'COMMON.LOCAL'
 991       include 'COMMON.NAMES'
 992       include 'COMMON.CHAIN'
 993       include 'COMMON.FFIELD'
 994       include 'COMMON.SBRIDGE'
 995       include 'COMMON.HEADER'
 996       include 'COMMON.CONTROL'
 997 c      include 'COMMON.DBASE'
 998 c      include 'COMMON.THREAD'
 999       include 'COMMON.TIME1'
1000 C Set up variable list.
1001       ntheta=nres-2
1002       nphi=nres-3
1003       nvar=ntheta+nphi
1004       nside=0
1005       do i=2,nres-1
1006         if (itype(i).ne.10) then
1007           nside=nside+1
1008           ialph(i,1)=nvar+nside
1009           ialph(nside,2)=i
1010         endif
1011       enddo
1012       if (indphi.gt.0) then
1013         nvar=nphi
1014       else if (indback.gt.0) then
1015         nvar=nphi+ntheta
1016       else
1017         nvar=nvar+2*nside
1018       endif
1019 cd    write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1020       return
1021       end
1022 c----------------------------------------------------------------------------
1023       subroutine gen_dist_constr
1024 C Generate CA distance constraints.
1025       implicit real*8 (a-h,o-z)
1026       include 'DIMENSIONS'
1027       include 'COMMON.IOUNITS'
1028       include 'COMMON.GEO'
1029       include 'COMMON.VAR'
1030       include 'COMMON.INTERACT'
1031       include 'COMMON.LOCAL'
1032       include 'COMMON.NAMES'
1033       include 'COMMON.CHAIN'
1034       include 'COMMON.FFIELD'
1035       include 'COMMON.SBRIDGE'
1036       include 'COMMON.HEADER'
1037       include 'COMMON.CONTROL'
1038 c      include 'COMMON.DBASE'
1039 c      include 'COMMON.THREAD'
1040       include 'COMMON.TIME1'
1041       dimension itype_pdb(maxres)
1042       common /pizda/ itype_pdb
1043       character*2 iden
1044 cd      print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1045 cd      write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1046 cd     & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1047 cd     & ' nsup',nsup
1048       do i=nstart_sup,nstart_sup+nsup-1
1049 cd      write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1050 cd     &    ' seq_pdb', restyp(itype_pdb(i))
1051         do j=i+2,nstart_sup+nsup-1
1052           nhpb=nhpb+1
1053           ihpb(nhpb)=i+nstart_seq-nstart_sup
1054           jhpb(nhpb)=j+nstart_seq-nstart_sup
1055           forcon(nhpb)=weidis
1056           dhpb(nhpb)=dist(i,j)
1057         enddo
1058       enddo
1059 cd      write (iout,'(a)') 'Distance constraints:'
1060 cd      do i=nss+1,nhpb
1061 cd        ii=ihpb(i)
1062 cd        jj=jhpb(i)
1063 cd        iden='CA'
1064 cd        if (ii.gt.nres) then
1065 cd          iden='SC'
1066 cd          ii=ii-nres
1067 cd          jj=jj-nres
1068 cd        endif
1069 cd        write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1070 cd     &  restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1071 cd     &  dhpb(i),forcon(i)
1072 cd      enddo
1073       return
1074       end
1075 c----------------------------------------------------------------------------
1076       subroutine csaread
1077       implicit real*8 (a-h,o-z)
1078       include 'DIMENSIONS'
1079       include 'COMMON.IOUNITS'
1080       include 'COMMON.GEO'
1081       include 'COMMON.CSA'
1082       include 'COMMON.BANK'
1083       include 'COMMON.CONTROL'
1084       character*80 ucase
1085       character*620 mcmcard
1086       call card_concat(mcmcard)
1087
1088       call readi(mcmcard,'NCONF',nconf,50)
1089       call readi(mcmcard,'NADD',nadd,0)
1090       call readi(mcmcard,'JSTART',jstart,1)
1091       call readi(mcmcard,'JEND',jend,1)
1092       call readi(mcmcard,'NSTMAX',nstmax,500000)
1093       call readi(mcmcard,'N0',n0,1)
1094       call readi(mcmcard,'N1',n1,6)
1095       call readi(mcmcard,'N2',n2,4)
1096       call readi(mcmcard,'N3',n3,0)
1097       call readi(mcmcard,'N4',n4,0)
1098       call readi(mcmcard,'N5',n5,0)
1099       call readi(mcmcard,'N6',n6,10)
1100       call readi(mcmcard,'N7',n7,0)
1101       call readi(mcmcard,'N8',n8,0)
1102       call readi(mcmcard,'N9',n9,0)
1103       call readi(mcmcard,'N14',n14,0)
1104       call readi(mcmcard,'N15',n15,0)
1105       call readi(mcmcard,'N16',n16,0)
1106       call readi(mcmcard,'N17',n17,0)
1107       call readi(mcmcard,'N18',n18,0)
1108
1109       vdisulf=(index(mcmcard,'DYNSS').gt.0)
1110
1111       call readi(mcmcard,'NDIFF',ndiff,2)
1112       call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1113       call readi(mcmcard,'IS1',is1,1)
1114       call readi(mcmcard,'IS2',is2,8)
1115       call readi(mcmcard,'NRAN0',nran0,4)
1116       call readi(mcmcard,'NRAN1',nran1,2)
1117       call readi(mcmcard,'IRR',irr,1)
1118       call readi(mcmcard,'NSEED',nseed,20)
1119       call readi(mcmcard,'NTOTAL',ntotal,10000)
1120       call reada(mcmcard,'CUT1',cut1,2.0d0)
1121       call reada(mcmcard,'CUT2',cut2,5.0d0)
1122       call reada(mcmcard,'ESTOP',estop,-300000.0d0)
1123       call readi(mcmcard,'ICMAX',icmax,1)
1124       call readi(mcmcard,'NBANKM',nbankm,400)
1125       call readi(mcmcard,'IUCUT',iucut,2)
1126       call readi(mcmcard,'IRESTART',irestart,0)
1127 c!bankt      call readi(mcmcard,'NBANKTM',ntbankm,0)
1128       ntbankm=0
1129 c!bankt
1130       call reada(mcmcard,'DELE',dele,20.0d0)
1131       call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1132       call readi(mcmcard,'IREF',iref,0)
1133       call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1134       call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1135       call readi(mcmcard,'NCONF_IN',nconf_in,0)
1136       call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1137       write (iout,*) "NCONF_IN",nconf_in
1138       tm_score=(index(mcmcard,'TMSCORE').gt.0)
1139       if (tm_score) write (iout,*) "Using TM_Score instead of DIFF",
1140      &      " for torsional angles"
1141       return
1142       end
1143
1144       subroutine read_minim
1145       implicit real*8 (a-h,o-z)
1146       include 'DIMENSIONS'
1147       include 'COMMON.MINIM'
1148       include 'COMMON.IOUNITS'
1149       character*80 ucase
1150       character*320 minimcard
1151       call card_concat(minimcard)
1152       call readi(minimcard,'MAXMIN',maxmin,2000)
1153       call readi(minimcard,'MAXFUN',maxfun,5000)
1154       call readi(minimcard,'MINMIN',minmin,maxmin)
1155       call readi(minimcard,'MINFUN',minfun,maxmin)
1156       call reada(minimcard,'TOLF',tolf,1.0D-2)
1157       call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1158       print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1159       print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1160       print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1161       write (iout,'(/80(1h*)/20x,a/80(1h*))')
1162      &         'Options in energy minimization:'
1163       write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1164      & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1165      & 'MinMin:',MinMin,' MinFun:',MinFun,
1166      & ' TolF:',TolF,' RTolF:',RTolF
1167       return
1168       end
1169 c----------------------------------------------------------------------------
1170       subroutine read_angles(kanal,*)
1171       implicit real*8 (a-h,o-z)
1172       include 'DIMENSIONS'
1173       include 'COMMON.GEO'
1174       include 'COMMON.VAR'
1175       include 'COMMON.CHAIN'
1176       include 'COMMON.IOUNITS'
1177       include 'COMMON.CONTROL'
1178 c Read angles from input
1179 c
1180        read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1181        read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1182        read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1183        read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1184
1185        do i=1,nres
1186 c 9/7/01 avoid 180 deg valence angle
1187         if (theta(i).gt.179.99d0) theta(i)=179.99d0
1188 c
1189         theta(i)=deg2rad*theta(i)
1190         phi(i)=deg2rad*phi(i)
1191         alph(i)=deg2rad*alph(i)
1192         omeg(i)=deg2rad*omeg(i)
1193        enddo
1194       return
1195    10 return1
1196       end
1197 c----------------------------------------------------------------------------
1198       subroutine reada(rekord,lancuch,wartosc,default)
1199       implicit none
1200       character*(*) rekord,lancuch
1201       double precision wartosc,default
1202       integer ilen,iread
1203       external ilen
1204       iread=index(rekord,lancuch)
1205       if (iread.eq.0) then
1206         wartosc=default
1207         return
1208       endif
1209       iread=iread+ilen(lancuch)+1
1210       read (rekord(iread:),*,err=10,end=10) wartosc
1211       return
1212   10  wartosc=default
1213       return
1214       end
1215 c----------------------------------------------------------------------------
1216       subroutine readi(rekord,lancuch,wartosc,default)
1217       implicit none
1218       character*(*) rekord,lancuch
1219       integer wartosc,default
1220       integer ilen,iread
1221       external ilen
1222       iread=index(rekord,lancuch)
1223       if (iread.eq.0) then
1224         wartosc=default
1225         return
1226       endif
1227       iread=iread+ilen(lancuch)+1
1228       read (rekord(iread:),*,err=10,end=10) wartosc
1229       return
1230   10  wartosc=default
1231       return
1232       end
1233 c----------------------------------------------------------------------------
1234       subroutine multreadi(rekord,lancuch,tablica,dim,default)
1235       implicit none
1236       integer dim,i
1237       integer tablica(dim),default
1238       character*(*) rekord,lancuch
1239       character*80 aux
1240       integer ilen,iread
1241       external ilen
1242       do i=1,dim
1243         tablica(i)=default
1244       enddo
1245       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1246       if (iread.eq.0) return
1247       iread=iread+ilen(lancuch)+1
1248       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1249    10 return
1250       end
1251 c----------------------------------------------------------------------------
1252       subroutine multreada(rekord,lancuch,tablica,dim,default)
1253       implicit none
1254       integer dim,i
1255       double precision tablica(dim),default
1256       character*(*) rekord,lancuch
1257       character*80 aux
1258       integer ilen,iread
1259       external ilen
1260       do i=1,dim
1261         tablica(i)=default
1262       enddo
1263       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1264       if (iread.eq.0) return
1265       iread=iread+ilen(lancuch)+1
1266       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1267    10 return
1268       end
1269 c----------------------------------------------------------------------------
1270       subroutine openunits
1271       implicit real*8 (a-h,o-z)
1272       include 'DIMENSIONS'
1273 #ifdef MPI
1274       include 'mpif.h'
1275       character*16 form,nodename
1276       integer nodelen
1277 #endif
1278       include 'COMMON.SETUP'
1279       include 'COMMON.IOUNITS'
1280 c      include 'COMMON.MD'
1281       include 'COMMON.CONTROL'
1282       integer lenpre,lenpot,ilen,lentmp
1283       external ilen
1284       character*3 out1file_text,ucase
1285       character*3 ll
1286       external ucase
1287 c      print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1288       call getenv_loc("PREFIX",prefix)
1289       pref_orig = prefix
1290       call getenv_loc("POT",pot)
1291       call getenv_loc("DIRTMP",tmpdir)
1292       call getenv_loc("CURDIR",curdir)
1293       call getenv_loc("OUT1FILE",out1file_text)
1294 c      print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1295       out1file_text=ucase(out1file_text)
1296       if (out1file_text(1:1).eq."Y") then
1297         out1file=.true.
1298       else
1299         out1file=fg_rank.gt.0
1300       endif
1301       lenpre=ilen(prefix)
1302       lenpot=ilen(pot)
1303       lentmp=ilen(tmpdir)
1304       if (lentmp.gt.0) then
1305           write (*,'(80(1h!))')
1306           write (*,'(a,19x,a,19x,a)') "!","  A T T E N T I O N  ","!"
1307           write (*,'(80(1h!))')
1308           write (*,*)"All output files will be on node /tmp directory."
1309 #ifdef MPI
1310         call  MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1311         if (me.eq.king) then
1312           write (*,*) "The master node is ",nodename
1313         else if (fg_rank.eq.0) then
1314           write (*,*) "I am the CG slave node ",nodename
1315         else
1316           write (*,*) "I am the FG slave node ",nodename
1317         endif
1318 #endif
1319         PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1320         lenpre = lentmp+lenpre+1
1321       endif
1322       entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
1323 C Get the names and open the input files
1324 #if defined(WINIFL) || defined(WINPGI)
1325       open(1,file=pref_orig(:ilen(pref_orig))//
1326      &  '.inp',status='old',readonly,shared)
1327        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1328 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1329 C Get parameter filenames and open the parameter files.
1330       call getenv_loc('BONDPAR',bondname)
1331       open (ibond,file=bondname,status='old',readonly,shared)
1332       call getenv_loc('THETPAR',thetname)
1333       open (ithep,file=thetname,status='old',readonly,shared)
1334 #ifndef CRYST_THETA
1335       call getenv_loc('THETPARPDB',thetname_pdb)
1336       open (ithep_pdb,file=thetname_pdb,status='old',readonly,shared)
1337 #endif
1338       call getenv_loc('ROTPAR',rotname)
1339       open (irotam,file=rotname,status='old',readonly,shared)
1340 #ifndef CRYST_SC
1341       call getenv_loc('ROTPARPDB',rotname_pdb)
1342       open (irotam_pdb,file=rotname_pdb,status='old',readonly,shared)
1343 #endif
1344       call getenv_loc('TORPAR',torname)
1345       open (itorp,file=torname,status='old',readonly,shared)
1346       call getenv_loc('TORDPAR',tordname)
1347       open (itordp,file=tordname,status='old',readonly,shared)
1348       call getenv_loc('FOURIER',fouriername)
1349       open (ifourier,file=fouriername,status='old',readonly,shared)
1350       call getenv_loc('ELEPAR',elename)
1351       open (ielep,file=elename,status='old',readonly,shared)
1352       call getenv_loc('SIDEPAR',sidename)
1353       open (isidep,file=sidename,status='old',readonly,shared)
1354 #elif (defined CRAY) || (defined AIX)
1355       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1356      &  action='read')
1357 c      print *,"Processor",myrank," opened file 1"
1358       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1359 c      print *,"Processor",myrank," opened file 9"
1360 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1361 C Get parameter filenames and open the parameter files.
1362       call getenv_loc('BONDPAR',bondname)
1363       open (ibond,file=bondname,status='old',action='read')
1364 c      print *,"Processor",myrank," opened file IBOND"
1365       call getenv_loc('THETPAR',thetname)
1366       open (ithep,file=thetname,status='old',action='read')
1367 c      print *,"Processor",myrank," opened file ITHEP"
1368 #ifndef CRYST_THETA
1369       call getenv_loc('THETPARPDB',thetname_pdb)
1370       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
1371 #endif
1372       call getenv_loc('ROTPAR',rotname)
1373       open (irotam,file=rotname,status='old',action='read')
1374 c      print *,"Processor",myrank," opened file IROTAM"
1375 #ifndef CRYST_SC
1376       call getenv_loc('ROTPARPDB',rotname_pdb)
1377       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
1378 #endif
1379       call getenv_loc('TORPAR',torname)
1380       open (itorp,file=torname,status='old',action='read')
1381 c      print *,"Processor",myrank," opened file ITORP"
1382       call getenv_loc('TORDPAR',tordname)
1383       open (itordp,file=tordname,status='old',action='read')
1384 c      print *,"Processor",myrank," opened file ITORDP"
1385       call getenv_loc('SCCORPAR',sccorname)
1386       open (isccor,file=sccorname,status='old',action='read')
1387 c      print *,"Processor",myrank," opened file ISCCOR"
1388       call getenv_loc('FOURIER',fouriername)
1389       open (ifourier,file=fouriername,status='old',action='read')
1390 c      print *,"Processor",myrank," opened file IFOURIER"
1391       call getenv_loc('ELEPAR',elename)
1392       open (ielep,file=elename,status='old',action='read')
1393 c      print *,"Processor",myrank," opened file IELEP"
1394       call getenv_loc('SIDEPAR',sidename)
1395       open (isidep,file=sidename,status='old',action='read')
1396 c      print *,"Processor",myrank," opened file ISIDEP"
1397 c      print *,"Processor",myrank," opened parameter files"
1398 #elif (defined G77)
1399       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
1400       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1401 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1402 C Get parameter filenames and open the parameter files.
1403       call getenv_loc('BONDPAR',bondname)
1404       open (ibond,file=bondname,status='old')
1405       call getenv_loc('THETPAR',thetname)
1406       open (ithep,file=thetname,status='old')
1407 #ifndef CRYST_THETA
1408       call getenv_loc('THETPARPDB',thetname_pdb)
1409       open (ithep_pdb,file=thetname_pdb,status='old')
1410 #endif
1411       call getenv_loc('ROTPAR',rotname)
1412       open (irotam,file=rotname,status='old')
1413 #ifndef CRYST_SC
1414       call getenv_loc('ROTPARPDB',rotname_pdb)
1415       open (irotam_pdb,file=rotname_pdb,status='old')
1416 #endif
1417       call getenv_loc('TORPAR',torname)
1418       open (itorp,file=torname,status='old')
1419       call getenv_loc('TORDPAR',tordname)
1420       open (itordp,file=tordname,status='old')
1421       call getenv_loc('SCCORPAR',sccorname)
1422       open (isccor,file=sccorname,status='old')
1423       call getenv_loc('FOURIER',fouriername)
1424       open (ifourier,file=fouriername,status='old')
1425       call getenv_loc('ELEPAR',elename)
1426       open (ielep,file=elename,status='old')
1427       call getenv_loc('SIDEPAR',sidename)
1428       open (isidep,file=sidename,status='old')
1429 #else
1430       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1431      &  readonly)
1432        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1433 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1434 C Get parameter filenames and open the parameter files.
1435       call getenv_loc('BONDPAR',bondname)
1436       open (ibond,file=bondname,status='old',readonly)
1437       call getenv_loc('THETPAR',thetname)
1438       open (ithep,file=thetname,status='old',readonly)
1439 #ifndef CRYST_THETA
1440       call getenv_loc('THETPARPDB',thetname_pdb)
1441       print *,"thetname_pdb ",thetname_pdb
1442       open (ithep_pdb,file=thetname_pdb,status='old',readonly)
1443       print *,ithep_pdb," opened"
1444 #endif
1445       call getenv_loc('ROTPAR',rotname)
1446       open (irotam,file=rotname,status='old',readonly)
1447 #ifndef CRYST_SC
1448       call getenv_loc('ROTPARPDB',rotname_pdb)
1449       open (irotam_pdb,file=rotname_pdb,status='old',readonly)
1450 #endif
1451       call getenv_loc('TORPAR',torname)
1452       open (itorp,file=torname,status='old',readonly)
1453       call getenv_loc('TORDPAR',tordname)
1454       open (itordp,file=tordname,status='old',readonly)
1455       call getenv_loc('SCCORPAR',sccorname)
1456       open (isccor,file=sccorname,status='old',readonly)
1457       call getenv_loc('FOURIER',fouriername)
1458       open (ifourier,file=fouriername,status='old',readonly)
1459       call getenv_loc('ELEPAR',elename)
1460       open (ielep,file=elename,status='old',readonly)
1461       call getenv_loc('SIDEPAR',sidename)
1462       open (isidep,file=sidename,status='old',readonly)
1463 #endif
1464 #ifndef OLDSCP
1465 C
1466 C 8/9/01 In the newest version SCp interaction constants are read from a file
1467 C Use -DOLDSCP to use hard-coded constants instead.
1468 C
1469       call getenv_loc('SCPPAR',scpname)
1470 #if defined(WINIFL) || defined(WINPGI)
1471       open (iscpp,file=scpname,status='old',readonly,shared)
1472 #elif (defined CRAY)  || (defined AIX)
1473       open (iscpp,file=scpname,status='old',action='read')
1474 #elif (defined G77)
1475       open (iscpp,file=scpname,status='old')
1476 #else
1477       open (iscpp,file=scpname,status='old',readonly)
1478 #endif
1479 #endif
1480       call getenv_loc('PATTERN',patname)
1481 #if defined(WINIFL) || defined(WINPGI)
1482       open (icbase,file=patname,status='old',readonly,shared)
1483 #elif (defined CRAY)  || (defined AIX)
1484       open (icbase,file=patname,status='old',action='read')
1485 #elif (defined G77)
1486       open (icbase,file=patname,status='old')
1487 #else
1488       open (icbase,file=patname,status='old',readonly)
1489 #endif
1490 #ifdef MPI
1491 C Open output file only for CG processes
1492 c      print *,"Processor",myrank," fg_rank",fg_rank
1493       if (fg_rank.eq.0) then
1494
1495       if (nodes.eq.1) then
1496         npos=3
1497       else
1498         npos = dlog10(dfloat(nodes-1))+1
1499       endif
1500       if (npos.lt.3) npos=3
1501       write (liczba,'(i1)') npos
1502       form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
1503      &  //')'
1504       write (liczba,form) me
1505       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
1506      &  liczba(:ilen(liczba))
1507       intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
1508      &  //'.int'
1509       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
1510      &  //'.pdb'
1511       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
1512      &  liczba(:ilen(liczba))//'.mol2'
1513       statname=prefix(:lenpre)//'_'//pot(:lenpot)//
1514      &  liczba(:ilen(liczba))//'.stat'
1515       if (lentmp.gt.0)
1516      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
1517      &      //liczba(:ilen(liczba))//'.stat')
1518       rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
1519      &  //'.rst'
1520 c      if(usampl) then
1521 c          qname=prefix(:lenpre)//'_'//pot(:lenpot)//
1522 c     & liczba(:ilen(liczba))//'.const'
1523 c      endif
1524
1525       endif
1526 #else
1527       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
1528       intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
1529       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
1530       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
1531       statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
1532       if (lentmp.gt.0)
1533      & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)//
1534      &    '.stat')
1535       rest2name=prefix(:ilen(prefix))//'.rst'
1536 c      if(usampl) then
1537 c         qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
1538 c      endif
1539 #endif
1540 #if defined(AIX) || defined(PGI)
1541       if (me.eq.king .or. .not. out1file)
1542      &   open(iout,file=outname,status='unknown')
1543 #ifdef DEBUG
1544       if (fg_rank.gt.0) then
1545         write (liczba,'(i3.3)') myrank/nfgtasks
1546         write (ll,'(bz,i3.3)') fg_rank
1547         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
1548      &   status='unknown')
1549       endif
1550 #endif
1551       if(me.eq.king) then
1552        open(igeom,file=intname,status='unknown',position='append')
1553        open(ipdb,file=pdbname,status='unknown')
1554        open(imol2,file=mol2name,status='unknown')
1555        open(istat,file=statname,status='unknown',position='append')
1556       else
1557 c1out       open(iout,file=outname,status='unknown')
1558       endif
1559 #else
1560       if (me.eq.king .or. .not.out1file)
1561      &    open(iout,file=outname,status='unknown')
1562 #ifdef DEBUG
1563       if (fg_rank.gt.0) then
1564         write (liczba,'(i3.3)') myrank/nfgtasks
1565         write (ll,'(bz,i3.3)') fg_rank
1566         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
1567      &   status='unknown')
1568       endif
1569 #endif
1570       if(me.eq.king) then
1571        open(igeom,file=intname,status='unknown',access='append')
1572        open(ipdb,file=pdbname,status='unknown')
1573        open(imol2,file=mol2name,status='unknown')
1574        open(istat,file=statname,status='unknown',access='append')
1575       else
1576 c1out       open(iout,file=outname,status='unknown')
1577       endif
1578 #endif
1579       csa_rbank=prefix(:lenpre)//'.CSA.rbank'
1580       csa_seed=prefix(:lenpre)//'.CSA.seed'
1581       csa_history=prefix(:lenpre)//'.CSA.history'
1582       csa_bank=prefix(:lenpre)//'.CSA.bank'
1583       csa_bank1=prefix(:lenpre)//'.CSA.bank1'
1584       csa_alpha=prefix(:lenpre)//'.CSA.alpha'
1585       csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
1586 c!bankt      csa_bankt=prefix(:lenpre)//'.CSA.bankt'
1587       csa_int=prefix(:lenpre)//'.int'
1588       csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
1589       csa_native_int=prefix(:lenpre)//'.CSA.native.int'
1590       csa_in=prefix(:lenpre)//'.CSA.in'
1591 c      print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
1592 C Write file names
1593       if (me.eq.king)then
1594       write (iout,'(80(1h-))')
1595       write (iout,'(30x,a)') "FILE ASSIGNMENT"
1596       write (iout,'(80(1h-))')
1597       write (iout,*) "Input file                      : ",
1598      &  pref_orig(:ilen(pref_orig))//'.inp'
1599       write (iout,*) "Output file                     : ",
1600      &  outname(:ilen(outname))
1601       write (iout,*)
1602       write (iout,*) "Sidechain potential file        : ",
1603      &  sidename(:ilen(sidename))
1604 #ifndef OLDSCP
1605       write (iout,*) "SCp potential file              : ",
1606      &  scpname(:ilen(scpname))
1607 #endif
1608       write (iout,*) "Electrostatic potential file    : ",
1609      &  elename(:ilen(elename))
1610       write (iout,*) "Cumulant coefficient file       : ",
1611      &  fouriername(:ilen(fouriername))
1612       write (iout,*) "Torsional parameter file        : ",
1613      &  torname(:ilen(torname))
1614       write (iout,*) "Double torsional parameter file : ",
1615      &  tordname(:ilen(tordname))
1616       write (iout,*) "SCCOR parameter file : ",
1617      &  sccorname(:ilen(sccorname))
1618       write (iout,*) "Bond & inertia constant file    : ",
1619      &  bondname(:ilen(bondname))
1620       write (iout,*) "Bending parameter file          : ",
1621      &  thetname(:ilen(thetname))
1622       write (iout,*) "Rotamer parameter file          : ",
1623      &  rotname(:ilen(rotname))
1624       write (iout,*) "Threading database              : ",
1625      &  patname(:ilen(patname))
1626       if (lentmp.ne.0)
1627      &write (iout,*)" DIRTMP                          : ",
1628      &  tmpdir(:lentmp)
1629       write (iout,'(80(1h-))')
1630       endif
1631       return
1632       end
1633 c----------------------------------------------------------------------------
1634       subroutine card_concat(card)
1635       implicit real*8 (a-h,o-z)
1636       include 'DIMENSIONS'
1637       include 'COMMON.IOUNITS'
1638       character*(*) card
1639       character*80 karta,ucase
1640       external ilen
1641       read (inp,'(a)') karta
1642       karta=ucase(karta)
1643       card=' '
1644       do while (karta(80:80).eq.'&')
1645         card=card(:ilen(card)+1)//karta(:79)
1646         read (inp,'(a)') karta
1647         karta=ucase(karta)
1648       enddo
1649       card=card(:ilen(card)+1)//karta
1650       return
1651       end
1652
1653       subroutine read_dist_constr
1654       implicit real*8 (a-h,o-z)
1655       include 'DIMENSIONS'
1656 #ifdef MPI
1657       include 'mpif.h'
1658 #endif
1659       include 'COMMON.SETUP'
1660       include 'COMMON.CONTROL'
1661       include 'COMMON.CHAIN'
1662       include 'COMMON.IOUNITS'
1663       include 'COMMON.SBRIDGE'
1664       integer ifrag_(2,100),ipair_(2,100)
1665       double precision wfrag_(100),wpair_(100)
1666       character*500 controlcard
1667       write (iout,*) "Calling read_dist_constr"
1668       write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
1669       call flush(iout)
1670       call card_concat(controlcard)
1671       call readi(controlcard,"NFRAG",nfrag_,0)
1672       call readi(controlcard,"NPAIR",npair_,0)
1673       call readi(controlcard,"NDIST",ndist_,0)
1674       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
1675       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
1676       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
1677       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
1678       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
1679 c      write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
1680 c      write (iout,*) "IFRAG"
1681 c      do i=1,nfrag_
1682 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
1683 c      enddo
1684 c      write (iout,*) "IPAIR"
1685 c      do i=1,npair_
1686 c        write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
1687 c      enddo
1688       call flush(iout)
1689       do i=1,nfrag_
1690         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
1691         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
1692      &    ifrag_(2,i)=nstart_sup+nsup-1
1693 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
1694         call flush(iout)
1695         if (wfrag_(i).gt.0.0d0) then
1696         do j=ifrag_(1,i),ifrag_(2,i)-1
1697           do k=j+1,ifrag_(2,i)
1698             write (iout,*) "j",j," k",k
1699             ddjk=dist(j,k)
1700             if (constr_dist.eq.1) then
1701             nhpb=nhpb+1
1702             ihpb(nhpb)=j
1703             jhpb(nhpb)=k
1704               dhpb(nhpb)=ddjk
1705             forcon(nhpb)=wfrag_(i)
1706             else if (constr_dist.eq.2) then
1707               if (ddjk.le.dist_cut) then
1708                 nhpb=nhpb+1
1709                 ihpb(nhpb)=j
1710                 jhpb(nhpb)=k
1711                 dhpb(nhpb)=ddjk
1712                 forcon(nhpb)=wfrag_(i)
1713               endif
1714             else
1715               nhpb=nhpb+1
1716               ihpb(nhpb)=j
1717               jhpb(nhpb)=k
1718               dhpb(nhpb)=ddjk
1719               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
1720             endif
1721 #ifdef MPI
1722             if (.not.out1file .or. me.eq.king)
1723      &      write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
1724      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1725 #else
1726             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
1727      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1728 #endif
1729           enddo
1730         enddo
1731         endif
1732       enddo
1733       do i=1,npair_
1734         if (wpair_(i).gt.0.0d0) then
1735         ii = ipair_(1,i)
1736         jj = ipair_(2,i)
1737         if (ii.gt.jj) then
1738           itemp=ii
1739           ii=jj
1740           jj=itemp
1741         endif
1742         do j=ifrag_(1,ii),ifrag_(2,ii)
1743           do k=ifrag_(1,jj),ifrag_(2,jj)
1744             nhpb=nhpb+1
1745             ihpb(nhpb)=j
1746             jhpb(nhpb)=k
1747             forcon(nhpb)=wpair_(i)
1748             dhpb(nhpb)=dist(j,k)
1749 #ifdef MPI
1750             if (.not.out1file .or. me.eq.king)
1751      &      write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1752      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1753 #else
1754             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1755      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1756 #endif
1757           enddo
1758         enddo
1759         endif
1760       enddo
1761       do i=1,ndist_
1762         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
1763         if (forcon(nhpb+1).gt.0.0d0) then
1764           nhpb=nhpb+1
1765           dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
1766 #ifdef MPI
1767           if (.not.out1file .or. me.eq.king)
1768      &    write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1769      &     nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1770 #else
1771           write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1772      &     nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1773 #endif
1774         endif
1775       enddo
1776       call flush(iout)
1777       return
1778       end
1779 c-------------------------------------------------------------------------------
1780 #ifdef WINIFL
1781       subroutine flush(iu)
1782       return
1783       end
1784 #endif
1785 #ifdef AIX
1786       subroutine flush(iu)
1787       call flush_(iu)
1788       return
1789       end
1790 #endif
1791 c------------------------------------------------------------------------------
1792       subroutine copy_to_tmp(source)
1793       include "DIMENSIONS"
1794       include "COMMON.IOUNITS"
1795       character*(*) source
1796       character* 256 tmpfile
1797       integer ilen
1798       external ilen
1799       logical ex
1800       tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
1801       inquire(file=tmpfile,exist=ex)
1802       if (ex) then
1803         write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
1804      &   " to temporary directory..."
1805         write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
1806         call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
1807       endif
1808       return
1809       end
1810 c------------------------------------------------------------------------------
1811       subroutine move_from_tmp(source)
1812       include "DIMENSIONS"
1813       include "COMMON.IOUNITS"
1814       character*(*) source
1815       integer ilen
1816       external ilen
1817       write (*,*) "Moving ",source(:ilen(source)),
1818      & " from temporary directory to working directory"
1819       write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
1820       call system("/bin/mv "//source(:ilen(source))//" "//curdir)
1821       return
1822       end
1823 c------------------------------------------------------------------------------
1824       subroutine random_init(seed)
1825 C
1826 C Initialize random number generator
1827 C
1828       implicit real*8 (a-h,o-z)
1829       include 'DIMENSIONS'
1830 #ifdef AMD64
1831       integer*8 iseedi8
1832 #endif
1833 #ifdef MPI
1834       include 'mpif.h'
1835       logical OKRandom, prng_restart
1836       real*8  r1
1837       integer iseed_array(4)
1838 #endif
1839       include 'COMMON.IOUNITS'
1840       include 'COMMON.TIME1'
1841 c      include 'COMMON.THREAD'
1842       include 'COMMON.SBRIDGE'
1843       include 'COMMON.CONTROL'
1844       include 'COMMON.MCM'
1845 c      include 'COMMON.MAP'
1846       include 'COMMON.HEADER'
1847 c      include 'COMMON.CSA'
1848       include 'COMMON.CHAIN'
1849 c      include 'COMMON.MUCA'
1850 c      include 'COMMON.MD'
1851       include 'COMMON.FFIELD'
1852       include 'COMMON.SETUP'
1853       iseed=-dint(dabs(seed))
1854       if (iseed.eq.0) then
1855         write (iout,'(/80(1h*)/20x,a/80(1h*))')
1856      &    'Random seed undefined. The program will stop.'
1857         write (*,'(/80(1h*)/20x,a/80(1h*))')
1858      &    'Random seed undefined. The program will stop.'
1859 #ifdef MPI
1860         call mpi_finalize(mpi_comm_world,ierr)
1861 #endif
1862         stop 'Bad random seed.'
1863       endif
1864 #ifdef MPI
1865       if (fg_rank.eq.0) then
1866       seed=seed*(me+1)+1
1867 #ifdef AMD64
1868       iseedi8=dint(seed)
1869       if(me.eq.king .or. .not. out1file)
1870      &  write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8
1871       write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8
1872       OKRandom = prng_restart(me,iseedi8)
1873 #else
1874       do i=1,4
1875        tmp=65536.0d0**(4-i)
1876        iseed_array(i) = dint(seed/tmp)
1877        seed=seed-iseed_array(i)*tmp
1878       enddo
1879       if(me.eq.king .or. .not. out1file)
1880      & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
1881      &                 (iseed_array(i),i=1,4)
1882       write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
1883      &                 (iseed_array(i),i=1,4)
1884       OKRandom = prng_restart(me,iseed_array)
1885 #endif
1886       if (OKRandom) then
1887         r1=ran_number(0.0D0,1.0D0)
1888         if(me.eq.king .or. .not. out1file)
1889      &   write (iout,*) 'ran_num',r1
1890         if (r1.lt.0.0d0) OKRandom=.false.
1891       endif
1892       if (.not.OKRandom) then
1893         write (iout,*) 'PRNG IS NOT WORKING!!!'
1894         print *,'PRNG IS NOT WORKING!!!'
1895         if (me.eq.0) then
1896          call flush(iout)
1897          call mpi_abort(mpi_comm_world,error_msg,ierr)
1898          stop
1899         else
1900          write (iout,*) 'too many processors for parallel prng'
1901          write (*,*) 'too many processors for parallel prng'
1902          call flush(iout)
1903          stop
1904         endif
1905       endif
1906       endif
1907 #else
1908       call vrndst(iseed)
1909       write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
1910 #endif
1911       return
1912       end