source/unres/src-HCD-5D/readrtns_CSA.F

   1       subroutine readrtns
   2       implicit none
   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       integer i,j
  13       logical file_exist
  14 C Read job setup parameters
  15       call read_control
  16 C Read force-field parameters except weights
  17       call parmread
  18 C Read control parameters for energy minimzation if required
  19       if (minim) call read_minim
  20 C Read MCM control parameters if required
  21       if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
  22 C Read MD control parameters if reqjuired
  23       if (modecalc.eq.12) call read_MDpar
  24 C Read MREMD control parameters if required
  25       if (modecalc.eq.14) then
  26          call read_MDpar
  27          call read_REMDpar
  28       endif
  29 C Read MUCA control parameters if required
  30       if (lmuca) call read_muca
  31 C Read CSA control parameters if required (from fort.40 if exists
  32 C otherwise from general input file)
  33       if (modecalc.eq.8) then
  34        inquire (file="fort.40",exist=file_exist)
  35        if (.not.file_exist) call csaread
  36       endif
  37 cfmc      if (modecalc.eq.10) call mcmfread
  38 c Read energy-term weights and disulfide parameters
  39       call weightread
  40 C Read molecule information, molecule geometry, energy-term weights, and
  41 C restraints if requested
  42       call molread
  43 C Print restraint information
  44 #ifdef MPI
  45       if (.not. out1file .or. me.eq.king) then
  46 #endif
  47       if (nhpb.gt.nss)  then
  48         write (iout,'(/a,i5,a/4a5,2a8,3a10,2a5)')
  49      &  "The following",nhpb-nss,
  50      &  " distance restraints have been imposed:",
  51      &  "   Nr"," res1"," res2"," beta","   d1","   d2","    k","    V",
  52      &  "  score"," type"
  53         do i=nss+1,nhpb
  54           write (iout,'(4i5,2f8.2,3f10.5,2i5)')i-nss,ihpb(i),jhpb(i),
  55      &  ibecarb(i),dhpb(i),dhpb1(i),forcon(i),fordepth(i),xlscore(i),
  56      &  irestr_type(i)
  57         enddo
  58       endif
  59       if (npeak.gt.0) then
  60         write (iout,'(/a,i5,a/4a5,2a8,3a10,2a5)')
  61      &  "The following",npeak,
  62      &  " NMR peak restraints have been imposed:",
  63      &  "   Nr"," res1"," res2"," beta","   d1","   d2","    k","    V",
  64      &  "  score"," type"," ipeak"
  65         do i=1,npeak
  66           do j=ipeak(1,i),ipeak(2,i)
  67             write (iout,'(5i5,2f8.2,2f10.5,i5)')i,j,ihpb_peak(j),
  68      &       jhpb_peak(j),ibecarb_peak(j),dhpb_peak(j),dhpb1_peak(j),
  69      &       forcon_peak(j),fordepth_peak(i),irestr_type_peak(j)
  70           enddo
  71         enddo
  72         write (iout,*) "The ipeak array"
  73         do i=1,npeak
  74           write (iout,'(3i5)' ) i,ipeak(1,i),ipeak(2,i)
  75         enddo
  76       endif
  77 #ifdef MPI
  78       endif
  79 #endif
  80 c      print *,"Processor",myrank," leaves READRTNS"
  81       return
  82       end
  83 C-------------------------------------------------------------------------------
  84       subroutine read_control
  85 C
  86 C Read contorl data
  87 C
  88       implicit none
  89       include 'DIMENSIONS'
  90 #ifdef MP
  91       include 'mpif.h'
  92       logical OKRandom, prng_restart
  93       double precision  r1
  94 #endif
  95       include 'COMMON.IOUNITS'
  96       include 'COMMON.TIME1'
  97       include 'COMMON.THREAD'
  98       include 'COMMON.SBRIDGE'
  99       include 'COMMON.CONTROL'
 100       include 'COMMON.SAXS'
 101       include 'COMMON.MCM'
 102       include 'COMMON.MAP'
 103       include 'COMMON.HEADER'
 104       include 'COMMON.CSA'
 105       include 'COMMON.CHAIN'
 106       include 'COMMON.MUCA'
 107       include 'COMMON.MD'
 108       include 'COMMON.FFIELD'
 109       include 'COMMON.INTERACT'
 110       include 'COMMON.SETUP'
 111       include 'COMMON.SPLITELE'
 112       include 'COMMON.SHIELD'
 113       include 'COMMON.GEO'
 114       integer i
 115       integer KDIAG,ICORFL,IXDR
 116       COMMON /MACHSW/ KDIAG,ICORFL,IXDR
 117       character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
 118       character*80 ucase
 119       character*320 controlcard
 120       double precision seed
 121
 122       nglob_csa=0
 123       eglob_csa=1d99
 124       nmin_csa=0
 125       read (INP,'(a)') titel
 126       call card_concat(controlcard)
 127 c      out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
 128 c      print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
 129       call reada(controlcard,'SEED',seed,0.0D0)
 130       call random_init(seed)
 131 C Set up the time limit (caution! The time must be input in minutes!)
 132       read_cart=index(controlcard,'READ_CART').gt.0
 133       out_cart=index(controlcard,'OUT_CART').gt.0
 134       out_int=index(controlcard,'OUT_INT').gt.0
 135       gmatout=index(controlcard,'GMATOUT').gt.0
 136       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
 137 C this variable with_theta_constr is the variable which allow to read and execute the
 138 C constrains on theta angles WITH_THETA_CONSTR is the keyword
 139       with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
 140       write (iout,*) "with_theta_constr ",with_theta_constr
 141       call readi(controlcard,'NSAXS',nsaxs,0)
 142       call readi(controlcard,'SAXS_MODE',saxs_mode,0)
 143       call reada(controlcard,'SCAL_RAD',scal_rad,1.0d0)
 144       call reada(controlcard,'SAXS_CUTOFF',saxs_cutoff,1.0d0)
 145       write (iout,*) "Number of SAXS restraints",NSAXS," SAXS_MODE",
 146      &   SAXS_MODE," SCAL_RAD",scal_rad,"SAXS_CUTOFF",saxs_cutoff
 147       call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
 148       read_homol_frag = index(controlcard,"READ_HOMOL_FRAG").gt.0
 149       out_template_coord = index(controlcard,"OUT_TEMPLATE_COORD").gt.0
 150       out_template_restr = index(controlcard,"OUT_TEMPLATE_RESTR").gt.0
 151       call readi(controlcard,'SYM',symetr,1)
 152       call reada(controlcard,'TIMLIM',timlim,2800.0D0) ! default 16 hours
 153       unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
 154       call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
 155 c      call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
 156 c      call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
 157 c      call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
 158 c      call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
 159 c      call reada(controlcard,'DRMS',drms,0.1D0)
 160 c      if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 161 c       write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
 162 c       write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
 163 c       write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
 164 c       write (iout,'(a,f10.1)')'DRMS    = ',drms
 165 cc       write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
 166 c       write (iout,'(a,f10.1)') 'Time limit (min):',timlim
 167 c      endif
 168       call readi(controlcard,'NZ_START',nz_start,0)
 169       call readi(controlcard,'NZ_END',nz_end,0)
 170 c      call readi(controlcard,'IZ_SC',iz_sc,0)
 171       timlim=60.0D0*timlim
 172       safety = 60.0d0*safety
 173       modecalc=0
 174       call readi(controlcard,"INTER_LIST_UPDATE",imatupdate,100)
 175       call reada(controlcard,"T_BATH",t_bath,300.0d0)
 176       minim=(index(controlcard,'MINIMIZE').gt.0)
 177       dccart=(index(controlcard,'CART').gt.0)
 178       overlapsc=(index(controlcard,'OVERLAP').gt.0)
 179       overlapsc=.not.overlapsc
 180       searchsc=(index(controlcard,'SEARCHSC').gt.0 .and.
 181      &  index(controlcard,'NOSEARCHSC').eq.0)
 182       sideadd=(index(controlcard,'SIDEADD').gt.0)
 183       energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
 184       mremd_dec=(index(controlcard,'MREMD_DEC').gt.0)
 185       outpdb=(index(controlcard,'PDBOUT').gt.0)
 186       outmol2=(index(controlcard,'MOL2OUT').gt.0)
 187       pdbref=(index(controlcard,'PDBREF').gt.0)
 188       refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
 189       indpdb=index(controlcard,'PDBSTART')
 190       AFMlog=(index(controlcard,'AFM'))
 191       selfguide=(index(controlcard,'SELFGUIDE'))
 192 c      print *,'AFMlog',AFMlog,selfguide,"KUPA"
 193       call readi(controlcard,'TUBEMOD',tubelog,0)
 194 c      write (iout,*) TUBElog,"TUBEMODE"
 195       call readi(controlcard,'IPRINT',iprint,0)
 196 C SHIELD keyword sets if the shielding effect of side-chains is used
 197 C 0 denots no shielding is used all peptide are equally despite the
 198 C solvent accesible area
 199 C 1 the newly introduced function
 200 C 2 reseved for further possible developement
 201       call readi(controlcard,'SHIELD',shield_mode,0)
 202 C      if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 203         write(iout,*) "shield_mode",shield_mode
 204 C      endif
 205       call readi(controlcard,'TORMODE',tor_mode,0)
 206 C      if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 207         write(iout,*) "torsional and valence angle mode",tor_mode
 208       call readi(controlcard,'MAXGEN',maxgen,10000)
 209       call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
 210       call readi(controlcard,"KDIAG",kdiag,0)
 211       call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
 212       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
 213      & write (iout,*) "RESCALE_MODE",rescale_mode
 214       split_ene=index(controlcard,'SPLIT_ENE').gt.0
 215       if (index(controlcard,'REGULAR').gt.0.0D0) then
 216         call reada(controlcard,'WEIDIS',weidis,0.1D0)
 217         modecalc=1
 218         refstr=.true.
 219       endif
 220       if (refstr) then
 221       if (index(controlcard,"CASC").gt.0) then
 222         iz_sc=1
 223 c      else if (index(controlcard,"CAONLY").gt.0) then
 224 c        iz_sc=0
 225       else if (index(controlcard,"SCONLY").gt.0) then
 226         iz_sc=2
 227       else
 228         iz_sc=0
 229 c        write (iout,*) "ERROR! Must specify CASC CAONLY or SCONLY when",
 230 c     &    " specifying REFSTR, PDBREF or REGULAR."
 231 c        stop
 232       endif
 233       endif
 234       if (index(controlcard,'CHECKGRAD').gt.0) then
 235         modecalc=5
 236         if (index(controlcard,' CART').gt.0) then
 237           icheckgrad=1
 238         elseif (index(controlcard,'CARINT').gt.0) then
 239           icheckgrad=2
 240         else
 241           icheckgrad=3
 242         endif
 243         call reada(controlcard,'DELTA',aincr,1.0d-7)
 244 c        write (iout,*) "icheckgrad",icheckgrad
 245       elseif (index(controlcard,'THREAD').gt.0) then
 246         modecalc=2
 247         call readi(controlcard,'THREAD',nthread,0)
 248         if (nthread.gt.0) then
 249           call reada(controlcard,'WEIDIS',weidis,0.1D0)
 250         else
 251           if (fg_rank.eq.0)
 252      &    write (iout,'(a)')'A number has to follow the THREAD keyword.'
 253           stop 'Error termination in Read_Control.'
 254         endif
 255       else if (index(controlcard,'MCMA').gt.0) then
 256         modecalc=3
 257       else if (index(controlcard,'MCEE').gt.0) then
 258         modecalc=6
 259       else if (index(controlcard,'MULTCONF').gt.0) then
 260         modecalc=4
 261       else if (index(controlcard,'MAP').gt.0) then
 262         modecalc=7
 263         call readi(controlcard,'MAP',nmap,0)
 264       else if (index(controlcard,'CSA').gt.0) then
 265         modecalc=8
 266 crc      else if (index(controlcard,'ZSCORE').gt.0) then
 267 crc
 268 crc  ZSCORE is rm from UNRES, modecalc=9 is available
 269 crc
 270 crc        modecalc=9
 271 cfcm      else if (index(controlcard,'MCMF').gt.0) then
 272 cfmc        modecalc=10
 273       else if (index(controlcard,'SOFTREG').gt.0) then
 274         modecalc=11
 275       else if (index(controlcard,'CHECK_BOND').gt.0) then
 276         modecalc=-1
 277       else if (index(controlcard,'TEST').gt.0) then
 278         modecalc=-2
 279       else if (index(controlcard,'MD').gt.0) then
 280         modecalc=12
 281       else if (index(controlcard,'RE ').gt.0) then
 282         modecalc=14
 283       endif
 284
 285       lmuca=index(controlcard,'MUCA').gt.0
 286       call readi(controlcard,'MUCADYN',mucadyn,0)
 287       call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
 288       if (lmuca .and. (me.eq.king .or. .not.out1file ))
 289      & then
 290        write (iout,*) 'MUCADYN=',mucadyn
 291        write (iout,*) 'MUCASMOOTH=',muca_smooth
 292       endif
 293
 294       iscode=index(controlcard,'ONE_LETTER')
 295       indphi=index(controlcard,'PHI')
 296       indback=index(controlcard,'BACK')
 297       iranconf=index(controlcard,'RAND_CONF')
 298       start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
 299       extconf=(index(controlcard,'EXTCONF').gt.0)
 300       if (start_from_model) then
 301         iranconf=0
 302         extconf=.false.
 303       endif
 304       i2ndstr=index(controlcard,'USE_SEC_PRED')
 305       gradout=index(controlcard,'GRADOUT').gt.0
 306       gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
 307 C  DISTCHAINMAX become obsolete for periodic boundry condition
 308       call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
 309 C Reading the dimensions of box in x,y,z coordinates
 310       call reada(controlcard,'BOXX',boxxsize,100.0d0)
 311       call reada(controlcard,'BOXY',boxysize,100.0d0)
 312       call reada(controlcard,'BOXZ',boxzsize,100.0d0)
 313       write(iout,*) "Periodic box dimensions",boxxsize,boxysize,boxzsize
 314 c Cutoff range for interactions
 315       call reada(controlcard,"R_CUT_INT",r_cut_int,25.0d0)
 316       call reada(controlcard,"R_CUT_RESPA",r_cut_respa,2.0d0)
 317       call reada(controlcard,"LAMBDA",rlamb,0.3d0)
 318       write (iout,*) "Cutoff on interactions",r_cut_int
 319       write (iout,*)
 320      & "Cutoff in switching short and long range interactions in RESPA",
 321      & r_cut_respa
 322       write (iout,*) "lambda in switch function",rlamb
 323       call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
 324       call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
 325       if (lipthick.gt.0.0d0) then
 326        bordliptop=(boxzsize+lipthick)/2.0
 327        bordlipbot=bordliptop-lipthick
 328 C      endif
 329       if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
 330      & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
 331       buflipbot=bordlipbot+lipbufthick
 332       bufliptop=bordliptop-lipbufthick
 333       if ((lipbufthick*2.0d0).gt.lipthick)
 334      &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
 335       endif
 336       write(iout,*) "bordliptop=",bordliptop
 337       write(iout,*) "bordlipbot=",bordlipbot
 338       write(iout,*) "bufliptop=",bufliptop
 339       write(iout,*) "buflipbot=",buflipbot
 340       write (iout,*) "SHIELD MODE",shield_mode
 341       if (TUBElog.gt.0) then
 342        call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
 343        call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
 344        call reada(controlcard,"RTUBE",tubeR0,0.0d0)
 345        call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
 346        call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
 347        call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
 348        buftubebot=bordtubebot+tubebufthick
 349        buftubetop=bordtubetop-tubebufthick
 350       endif
 351       if (me.eq.king .or. .not.out1file )
 352      &  write (iout,*) "DISTCHAINMAX",distchainmax
 353
 354       if(me.eq.king.or..not.out1file)
 355      & write (iout,'(2a)') diagmeth(kdiag),
 356      &  ' routine used to diagonalize matrices.'
 357       return
 358       end
 359 c--------------------------------------------------------------------------
 360       subroutine read_REMDpar
 361 C
 362 C Read REMD settings
 363 C
 364       implicit none
 365       include 'DIMENSIONS'
 366       include 'COMMON.IOUNITS'
 367       include 'COMMON.TIME1'
 368       include 'COMMON.MD'
 369 #ifndef LANG0
 370       include 'COMMON.LANGEVIN'
 371 #else
 372 #ifdef FIVEDIAG
 373       include 'COMMON.LANGEVIN.lang0.5diag'
 374 #else
 375       include 'COMMON.LANGEVIN.lang0'
 376 #endif
 377 #endif
 378       include 'COMMON.INTERACT'
 379       include 'COMMON.NAMES'
 380       include 'COMMON.GEO'
 381       include 'COMMON.REMD'
 382       include 'COMMON.CONTROL'
 383       include 'COMMON.SETUP'
 384       character*80 ucase
 385       character*320 controlcard
 386       character*3200 controlcard1
 387       integer iremd_m_total,i
 388
 389       if(me.eq.king.or..not.out1file)
 390      & write (iout,*) "REMD setup"
 391
 392       call card_concat(controlcard)
 393       call readi(controlcard,"NREP",nrep,3)
 394       call readi(controlcard,"NSTEX",nstex,1000)
 395       call reada(controlcard,"RETMIN",retmin,10.0d0)
 396       call reada(controlcard,"RETMAX",retmax,1000.0d0)
 397       mremdsync=(index(controlcard,'SYNC').gt.0)
 398       call readi(controlcard,"NSYN",i_sync_step,100)
 399       restart1file=(index(controlcard,'REST1FILE').gt.0)
 400       traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
 401       call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
 402       if(max_cache_traj_use.gt.max_cache_traj)
 403      &           max_cache_traj_use=max_cache_traj
 404       if(me.eq.king.or..not.out1file) then
 405 cd       if (traj1file) then
 406 crc caching is in testing - NTWX is not ignored
 407 cd        write (iout,*) "NTWX value is ignored"
 408 cd        write (iout,*) "  trajectory is stored to one file by master"
 409 cd        write (iout,*) "  before exchange at NSTEX intervals"
 410 cd       endif
 411        write (iout,*) "NREP= ",nrep
 412        write (iout,*) "NSTEX= ",nstex
 413        write (iout,*) "SYNC= ",mremdsync
 414        write (iout,*) "NSYN= ",i_sync_step
 415        write (iout,*) "TRAJCACHE= ",max_cache_traj_use
 416       endif
 417       remd_tlist=.false.
 418       if (index(controlcard,'TLIST').gt.0) then
 419          remd_tlist=.true.
 420          call card_concat(controlcard1)
 421          read(controlcard1,*) (remd_t(i),i=1,nrep)
 422          if(me.eq.king.or..not.out1file)
 423      &    write (iout,*)'tlist',(remd_t(i),i=1,nrep)
 424       endif
 425       remd_mlist=.false.
 426       if (index(controlcard,'MLIST').gt.0) then
 427          remd_mlist=.true.
 428          call card_concat(controlcard1)
 429          read(controlcard1,*) (remd_m(i),i=1,nrep)
 430          if(me.eq.king.or..not.out1file) then
 431           write (iout,*)'mlist',(remd_m(i),i=1,nrep)
 432           iremd_m_total=0
 433           do i=1,nrep
 434            iremd_m_total=iremd_m_total+remd_m(i)
 435           enddo
 436            write (iout,*) 'Total number of replicas ',iremd_m_total
 437           endif
 438          endif
 439         if (adaptive) then
 440           write (iout,*)
 441      &      "Adaptive (PMF-biased) umbrella sampling will be run"
 442           call PMFread
 443         endif
 444       if(me.eq.king.or..not.out1file)
 445      &   write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
 446       return
 447       end
 448 c--------------------------------------------------------------------------
 449       subroutine read_MDpar
 450 C
 451 C Read MD settings
 452 C
 453       implicit none
 454       include 'DIMENSIONS'
 455       include 'COMMON.IOUNITS'
 456       include 'COMMON.TIME1'
 457       include 'COMMON.MD'
 458       include 'COMMON.QRESTR'
 459 #ifndef LANG0
 460       include 'COMMON.LANGEVIN'
 461 #else
 462 #ifdef FIVEDIAG
 463       include 'COMMON.LANGEVIN.lang0.5diag'
 464 #else
 465       include 'COMMON.LANGEVIN.lang0'
 466 #endif
 467 #endif
 468       include 'COMMON.INTERACT'
 469       include 'COMMON.NAMES'
 470       include 'COMMON.GEO'
 471       include 'COMMON.SETUP'
 472       include 'COMMON.CONTROL'
 473       include 'COMMON.SPLITELE'
 474       include 'COMMON.FFIELD'
 475       character*80 ucase
 476       character*320 controlcard
 477       integer i
 478       double precision eta
 479
 480       call card_concat(controlcard)
 481       call readi(controlcard,"NSTEP",n_timestep,1000000)
 482       call readi(controlcard,"NTWE",ntwe,100)
 483       call readi(controlcard,"NTWX",ntwx,1000)
 484       call reada(controlcard,"DT",d_time,1.0d-1)
 485       call reada(controlcard,"DVMAX",dvmax,2.0d1)
 486       call reada(controlcard,"DAMAX",damax,1.0d1)
 487       call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
 488       call readi(controlcard,"LANG",lang,0)
 489       RESPA = index(controlcard,"RESPA") .gt. 0
 490       call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
 491       ntime_split0=ntime_split
 492       call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
 493       ntime_split0=ntime_split
 494 c      call reada(controlcard,"R_CUT",r_cut,2.0d0)
 495 c      call reada(controlcard,"LAMBDA",rlamb,0.3d0)
 496       rest = index(controlcard,"REST").gt.0
 497       tbf = index(controlcard,"TBF").gt.0
 498       usampl = index(controlcard,"USAMPL").gt.0
 499       scale_umb =  index(controlcard,"SCALE_UMB").gt.0
 500       adaptive = index(controlcard,"ADAPTIVE").gt.0
 501       mdpdb = index(controlcard,"MDPDB").gt.0
 502       call reada(controlcard,"T_BATH",t_bath,300.0d0)
 503       call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
 504       call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
 505       call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
 506       if (count_reset_moment.eq.0) count_reset_moment=1000000000
 507       call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
 508       reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
 509       reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
 510       if (count_reset_vel.eq.0) count_reset_vel=1000000000
 511       large = index(controlcard,"LARGE").gt.0
 512       print_compon = index(controlcard,"PRINT_COMPON").gt.0
 513       rattle = index(controlcard,"RATTLE").gt.0
 514       preminim = index(controlcard,"PREMINIM").gt.0
 515       if (iranconf.gt.0 .or. indpdb.gt.0 .or. start_from_model) then
 516         preminim=.true.
 517         dccart=.true.
 518       endif
 519       write (iout,*) "PREMINIM ",preminim
 520       if (preminim) then
 521         if (index(controlcard,'CART').gt.0) dccart=.true.
 522         write (iout,*) "dccart ",dccart
 523         write (iout,*) "read_minim ",index(controlcard,'READ_MINIM_PAR')
 524         if (index(controlcard,'READ_MINIM_PAR').gt.0) call read_minim
 525       endif
 526 c  if performing umbrella sampling, fragments constrained are read from the fragment file
 527       nset=0
 528       if(usampl) then
 529         write (iout,*) "Umbrella sampling will be run"
 530         if (scale_umb.and.adaptive) then
 531           write (iout,*) "ADAPTIVE and SCALE_UMB are mutually exclusive"
 532           write (iout,*) "Select one of those and re-run the job."
 533           stop
 534         endif
 535         if (scale_umb) write (iout,*)
 536      &"Umbrella-restraint force constants will be scaled by temperature"
 537         call read_fragments
 538       endif
 539
 540       if(me.eq.king.or..not.out1file) then
 541        write (iout,*)
 542        write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
 543        write (iout,*)
 544        write (iout,'(a)') "The units are:"
 545        write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
 546        write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
 547      &  " acceleration: angstrom/(48.9 fs)**2"
 548        write (iout,'(a)') "energy: kcal/mol, temperature: K"
 549        write (iout,*)
 550        write (iout,'(a60,i10)') "Number of time steps:",n_timestep
 551        write (iout,'(a60,f10.5,a)')
 552      &  "Initial time step of numerical integration:",d_time,
 553      &  " natural units"
 554        write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
 555        write (iout,'(a60,f10.5,a)') "Cutoff on interactions",r_cut_int,
 556      &   " A"
 557        write(iout,'(a60,i5)')"Frequency of updating interaction list",
 558      &   imatupdate
 559        if (RESPA) then
 560         write (iout,'(2a,i4,a)')
 561      &    "A-MTS algorithm used; initial time step for fast-varying",
 562      &    " short-range forces split into",ntime_split," steps."
 563         write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
 564      &   r_cut_respa," lambda",rlamb
 565        endif
 566        write (iout,'(2a,f10.5)')
 567      &  "Maximum acceleration threshold to reduce the time step",
 568      &  "/increase split number:",damax
 569        write (iout,'(2a,f10.5)')
 570      &  "Maximum predicted energy drift to reduce the timestep",
 571      &  "/increase split number:",edriftmax
 572        write (iout,'(a60,f10.5)')
 573      & "Maximum velocity threshold to reduce velocities:",dvmax
 574        write (iout,'(a60,i10)') "Frequency of property output:",ntwe
 575        write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
 576        if (rattle) write (iout,'(a60)')
 577      &  "Rattle algorithm used to constrain the virtual bonds"
 578       endif
 579       reset_fricmat=1000
 580       if (lang.gt.0) then
 581         call reada(controlcard,"ETAWAT",etawat,0.8904d0)
 582         call reada(controlcard,"RWAT",rwat,1.4d0)
 583         call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
 584         surfarea=index(controlcard,"SURFAREA").gt.0
 585         call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
 586         if(me.eq.king.or..not.out1file)then
 587          write (iout,'(/a,$)') "Langevin dynamics calculation"
 588          if (lang.eq.1) then
 589           write (iout,'(a/)')
 590      &      " with direct integration of Langevin equations"
 591          else if (lang.eq.2) then
 592           write (iout,'(a/)') " with TINKER stochasic MD integrator"
 593          else if (lang.eq.3) then
 594           write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
 595          else if (lang.eq.4) then
 596           write (iout,'(a/)') " in overdamped mode"
 597          else
 598           write (iout,'(//a,i5)')
 599      &      "=========== ERROR: Unknown Langevin dynamics mode:",lang
 600           stop
 601          endif
 602          write (iout,'(a60,f10.5)') "Temperature:",t_bath
 603          write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
 604          write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
 605          write (iout,'(a60,f10.5)')
 606      &   "Scaling factor of the friction forces:",scal_fric
 607          if (surfarea) write (iout,'(2a,i10,a)')
 608      &     "Friction coefficients will be scaled by solvent-accessible",
 609      &     " surface area every",reset_fricmat," steps."
 610         endif
 611 c Calculate friction coefficients and bounds of stochastic forces
 612         eta=6*pi*cPoise*etawat
 613         if(me.eq.king.or..not.out1file)
 614      &   write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
 615      &   ,eta
 616         gamp=scal_fric*(pstok+rwat)*eta
 617         stdfp=dsqrt(2*Rb*t_bath/d_time)
 618         do i=1,ntyp
 619           gamsc(i)=scal_fric*(restok(i)+rwat)*eta
 620           stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
 621         enddo
 622         if(me.eq.king.or..not.out1file)then
 623          write (iout,'(/2a/)')
 624      &   "Radii of site types and friction coefficients and std's of",
 625      &   " stochastic forces of fully exposed sites"
 626          write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
 627          do i=1,ntyp
 628           write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
 629      &     gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
 630          enddo
 631         endif
 632       else if (tbf) then
 633         if(me.eq.king.or..not.out1file)then
 634          write (iout,'(a)') "Berendsen bath calculation"
 635          write (iout,'(a60,f10.5)') "Temperature:",t_bath
 636          write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
 637          if (reset_moment)
 638      &   write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
 639      &   count_reset_moment," steps"
 640          if (reset_vel)
 641      &    write (iout,'(a,i10,a)')
 642      &    "Velocities will be reset at random every",count_reset_vel,
 643      &   " steps"
 644         endif
 645       else
 646         if(me.eq.king.or..not.out1file)
 647      &   write (iout,'(a31)') "Microcanonical mode calculation"
 648       endif
 649       if(me.eq.king.or..not.out1file)then
 650        if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
 651        if (usampl) then
 652           write(iout,*) "MD running with constraints."
 653           write(iout,*) "Equilibration time ", eq_time, " mtus."
 654           write(iout,*) "Constraining ", nfrag," fragments."
 655           write(iout,*) "Length of each fragment, weight and q0:"
 656           do iset=1,nset
 657            write (iout,*) "Set of restraints #",iset
 658            do i=1,nfrag
 659               write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
 660      &           ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
 661            enddo
 662            write(iout,*) "constraints between ", npair, "fragments."
 663            write(iout,*) "constraint pairs, weights and q0:"
 664            do i=1,npair
 665             write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
 666      &             ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
 667            enddo
 668            write(iout,*) "angle constraints within ", nfrag_back,
 669      &      "backbone fragments."
 670            if (loc_qlike) then
 671            write(iout,*) "fragment, weights, q0:"
 672            do i=1,nfrag_back
 673             write(iout,'(2i5,3(f8.1,f8.2))') ifrag_back(1,i,iset),
 674      &         ifrag_back(2,i,iset),
 675      &         wfrag_back(1,i,iset),qin_back(1,i,iset),
 676      &         wfrag_back(2,i,iset),qin_back(2,i,iset),
 677      &         wfrag_back(3,i,iset),qin_back(3,i,iset)
 678            enddo
 679            else
 680            write(iout,*) "fragment, weights:"
 681            do i=1,nfrag_back
 682             write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
 683      &         ifrag_back(2,i,iset),wfrag_back(1,i,iset),
 684      &         wfrag_back(2,i,iset),wfrag_back(3,i,iset)
 685            enddo
 686            endif
 687           enddo
 688         iset=mod(kolor,nset)+1
 689        endif
 690       endif
 691       if(me.eq.king.or..not.out1file)
 692      & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
 693       return
 694       end
 695 c------------------------------------------------------------------------------
 696       subroutine molread
 697 C
 698 C Read molecular data.
 699 C
 700       implicit none
 701       include 'DIMENSIONS'
 702 #ifdef MPI
 703       include 'mpif.h'
 704       integer error_msg,ierror,ierr,ierrcode
 705 #endif
 706       include 'COMMON.IOUNITS'
 707       include 'COMMON.GEO'
 708       include 'COMMON.VAR'
 709       include 'COMMON.INTERACT'
 710       include 'COMMON.LOCAL'
 711       include 'COMMON.NAMES'
 712       include 'COMMON.CHAIN'
 713       include 'COMMON.FFIELD'
 714       include 'COMMON.SBRIDGE'
 715       include 'COMMON.HEADER'
 716       include 'COMMON.CONTROL'
 717       include 'COMMON.SAXS'
 718       include 'COMMON.DBASE'
 719       include 'COMMON.THREAD'
 720       include 'COMMON.CONTACTS'
 721       include 'COMMON.TORCNSTR'
 722       include 'COMMON.TIME1'
 723       include 'COMMON.BOUNDS'
 724       include 'COMMON.MD'
 725       include 'COMMON.SETUP'
 726       include 'COMMON.SHIELD'
 727       character*4 sequence(maxres)
 728       integer rescode
 729       double precision x(maxvar)
 730       character*256 pdbfile
 731       character*400 weightcard
 732       character*80 weightcard_t,ucase
 733       integer itype_pdb(maxres)
 734       common /pizda/ itype_pdb
 735       logical seq_comp,fail
 736       double precision energia(0:n_ene)
 737       double precision secprob(3,maxdih_constr)
 738       double precision co
 739       double precision phihel,phibet,sigmahel,sigmabet
 740       integer iti,nsi,maxsi
 741       integer ilen
 742       external ilen
 743       integer iperm,tperm
 744       integer i,j,ii,k,l,itrial,itmp,i1,i2,it1,it2,nres_temp,itemp
 745       double precision sumv
 746 C
 747 C Read PDB structure if applicable
 748 C
 749       if (indpdb.gt.0 .or. pdbref) then
 750         read(inp,'(a)') pdbfile
 751         if(me.eq.king.or..not.out1file)
 752      &   write (iout,'(2a)') 'PDB data will be read from file ',
 753      &   pdbfile(:ilen(pdbfile))
 754         open(ipdbin,file=pdbfile,status='old',err=33)
 755         goto 34
 756   33    write (iout,'(a)') 'Error opening PDB file.'
 757         stop
 758   34    continue
 759         call readpdb
 760         do i=1,2*nres
 761           do j=1,3
 762             crefjlee(j,i)=c(j,i)
 763           enddo
 764         enddo
 765 #ifdef DEBUG
 766         do i=1,nres
 767           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
 768      &      (crefjlee(j,i+nres),j=1,3)
 769         enddo
 770 #endif
 771         if(me.eq.king.or..not.out1file)
 772      &   write (iout,'(a,i3,a,i3)')'nsup=',nsup,
 773      &   ' nstart_sup=',nstart_sup
 774         do i=1,nres
 775           itype_pdb(i)=itype(i)
 776         enddo
 777         close (ipdbin)
 778         nnt=nstart_sup
 779         nct=nstart_sup+nsup-1
 780         call contact(.false.,ncont_ref,icont_ref,co)
 781
 782         if (sideadd) then
 783          if(me.eq.king.or..not.out1file)
 784      &    write(iout,*)'Adding sidechains'
 785          maxsi=1000
 786          do i=2,nres-1
 787           iti=itype(i)
 788           if (iti.ne.10 .and. itype(i).ne.ntyp1) then
 789             nsi=0
 790             fail=.true.
 791             do while (fail.and.nsi.le.maxsi)
 792               call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
 793               nsi=nsi+1
 794             enddo
 795             if(fail) write(iout,*)'Adding sidechain failed for res ',
 796      &              i,' after ',nsi,' trials'
 797           endif
 798          enddo
 799         endif
 800       endif
 801       if (indpdb.eq.0) then
 802 C Read sequence if not taken from the pdb file.
 803         read (inp,*) nres
 804 c        print *,'nres=',nres
 805         if (iscode.gt.0) then
 806           read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
 807         else
 808           read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
 809         endif
 810 C Convert sequence to numeric code
 811         do i=1,nres
 812           itype(i)=rescode(i,sequence(i),iscode)
 813         enddo
 814       endif
 815 c      print *,nres
 816 c      print '(20i4)',(itype(i),i=1,nres)
 817       do i=1,nres
 818 #ifdef PROCOR
 819         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
 820 #else
 821         if (itype(i).eq.ntyp1) then
 822 #endif
 823           itel(i)=0
 824 #ifdef PROCOR
 825         else if (iabs(itype(i+1)).ne.20) then
 826 #else
 827         else if (iabs(itype(i)).ne.20) then
 828 #endif
 829           itel(i)=1
 830         else
 831           itel(i)=2
 832         endif
 833       enddo
 834       if(me.eq.king.or..not.out1file)then
 835        write (iout,*) "ITEL"
 836        do i=1,nres-1
 837          write (iout,*) i,itype(i),itel(i)
 838        enddo
 839 c       print *,'Call Read_Bridge.'
 840       endif
 841       nnt=1
 842       nct=nres
 843 cd      print *,'NNT=',NNT,' NCT=',NCT
 844       call seq2chains(nres,itype,nchain,chain_length,chain_border,
 845      &  ireschain)
 846       chain_border1(1,1)=1
 847       chain_border1(2,1)=chain_border(2,1)+1
 848       do i=2,nchain-1
 849         chain_border1(1,i)=chain_border(1,i)-1
 850         chain_border1(2,i)=chain_border(2,i)+1
 851       enddo
 852       if (nchain.gt.1) chain_border1(1,nchain)=chain_border(1,nchain)-1
 853       chain_border1(2,nchain)=nres
 854       write(iout,*) "nres",nres," nchain",nchain
 855       do i=1,nchain
 856         write(iout,*)"chain",i,chain_length(i),chain_border(1,i),
 857      &    chain_border(2,i),chain_border1(1,i),chain_border1(2,i)
 858       enddo
 859       call chain_symmetry(nchain,nres,itype,chain_border,
 860      &    chain_length,npermchain,tabpermchain)
 861 c      do i=1,nres
 862 c        write(iout,*) i,(tperm(ireschain(i),ii,tabpermchain),
 863 c     &    ii=1,npermchain)
 864 c      enddo
 865       write(iout,*) "residue permutations"
 866       do i=1,nres
 867         write(iout,*) i,(iperm(i,ii),ii=1,npermchain)
 868       enddo
 869       call flush(iout)
 870       if (itype(1).eq.ntyp1) nnt=2
 871       if (itype(nres).eq.ntyp1) nct=nct-1
 872       write (iout,*) "nnt",nnt," nct",nct
 873       call flush(iout)
 874 #ifdef DFA
 875       if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
 876      &            wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
 877        call init_dfa_vars
 878 c       print*, 'init_dfa_vars finished!'
 879        call read_dfa_info
 880 c       print*, 'read_dfa_info finished!'
 881       endif
 882 #endif
 883       if (pdbref) then
 884         if(me.eq.king.or..not.out1file)
 885      &   write (iout,'(a,i3)') 'nsup=',nsup
 886         nstart_seq=nnt
 887         if (nsup.le.(nct-nnt+1)) then
 888           do i=0,nct-nnt+1-nsup
 889             if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
 890               nstart_seq=nnt+i
 891               goto 111
 892             endif
 893           enddo
 894           write (iout,'(a)')
 895      &            'Error - sequences to be superposed do not match.'
 896           stop
 897         else
 898           do i=0,nsup-(nct-nnt+1)
 899             if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
 900      &      then
 901               nstart_sup=nstart_sup+i
 902               nsup=nct-nnt+1
 903               goto 111
 904             endif
 905           enddo
 906           write (iout,'(a)')
 907      &            'Error - sequences to be superposed do not match.'
 908         endif
 909   111   continue
 910         if (nsup.eq.0) nsup=nct-nnt
 911         if (nstart_sup.eq.0) nstart_sup=nnt
 912         if (nstart_seq.eq.0) nstart_seq=nnt
 913         if(me.eq.king.or..not.out1file)
 914      &   write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
 915      &                 ' nstart_seq=',nstart_seq
 916       endif
 917       call read_bridge
 918 C 8/13/98 Set limits to generating the dihedral angles
 919       do i=1,nres
 920         phibound(1,i)=-pi
 921         phibound(2,i)=pi
 922       enddo
 923       read (inp,*) ndih_constr
 924       write (iout,*) "ndih_constr",ndih_constr
 925       if (ndih_constr.gt.0) then
 926         raw_psipred=.false.
 927 C        read (inp,*) ftors
 928         read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
 929      &  i=1,ndih_constr)
 930 #ifdef MPI
 931         if(me.eq.king.or..not.out1file)then
 932 #endif
 933          write (iout,*)
 934      &   'There are',ndih_constr,' restraints on gamma angles.'
 935          do i=1,ndih_constr
 936           write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
 937      &    ftors(i)
 938          enddo
 939 #ifdef MPI
 940         endif
 941 #endif
 942         do i=1,ndih_constr
 943           phi0(i)=deg2rad*phi0(i)
 944           drange(i)=deg2rad*drange(i)
 945         enddo
 946 C        if(me.eq.king.or..not.out1file)
 947 C     &   write (iout,*) 'FTORS',ftors
 948         do i=1,ndih_constr
 949           ii = idih_constr(i)
 950           phibound(1,ii) = phi0(i)-drange(i)
 951           phibound(2,ii) = phi0(i)+drange(i)
 952         enddo
 953 #ifdef MPI
 954         if (me.eq.king .or. .not.out1file) then
 955 #endif
 956          write (iout,'(a)') 'Boundaries in gamma angle sampling:'
 957          do i=1,nres
 958            write (iout,'(a3,i5,2f10.1)')
 959      &    restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
 960          enddo
 961 #ifdef MP
 962         endif
 963 #endif
 964       else if (ndih_constr.lt.0) then
 965         raw_psipred=.true.
 966         call card_concat(weightcard)
 967         call reada(weightcard,"PHIHEL",phihel,50.0D0)
 968         call reada(weightcard,"PHIBET",phibet,180.0D0)
 969         call reada(weightcard,"SIGMAHEL",sigmahel,30.0d0)
 970         call reada(weightcard,"SIGMABET",sigmabet,40.0d0)
 971         call reada(weightcard,"WDIHC",wdihc,0.591D0)
 972         write (iout,*) "Weight of dihedral angle restraints",wdihc
 973         read(inp,'(9x,3f7.3)')
 974      &     (secprob(1,i),secprob(2,i),secprob(3,i),i=nnt,nct)
 975         write (iout,*) "The secprob array"
 976         do i=nnt,nct
 977           write (iout,'(i5,3f8.3)') i,(secprob(j,i),j=1,3)
 978         enddo
 979         ndih_constr=0
 980         do i=nnt+3,nct
 981           if (itype(i-3).ne.ntyp1 .and. itype(i-2).ne.ntyp1
 982      &    .and. itype(i-1).ne.ntyp1 .and. itype(i).ne.ntyp1) then
 983             ndih_constr=ndih_constr+1
 984             idih_constr(ndih_constr)=i
 985             sumv=0.0d0
 986             do j=1,3
 987               vpsipred(j,ndih_constr)=secprob(j,i-1)*secprob(j,i-2)
 988               sumv=sumv+vpsipred(j,ndih_constr)
 989             enddo
 990             do j=1,3
 991               vpsipred(j,ndih_constr)=vpsipred(j,ndih_constr)/sumv
 992             enddo
 993             phibound(1,ndih_constr)=phihel*deg2rad
 994             phibound(2,ndih_constr)=phibet*deg2rad
 995             sdihed(1,ndih_constr)=sigmahel*deg2rad
 996             sdihed(2,ndih_constr)=sigmabet*deg2rad
 997           endif
 998         enddo
 999 #ifdef MPI
1000         if(me.eq.king.or..not.out1file)then
1001 #endif
1002          write (iout,*)
1003      &   'There are',ndih_constr,
1004      &   ' bimodal restraints on gamma angles.'
1005          do i=1,ndih_constr
1006            write(iout,'(i5,1x,a4,i5,1h-,a4,i5,4f8.3,3f10.5)') i,
1007      &      restyp(itype(idih_constr(i)-2)),idih_constr(i)-2,
1008      &      restyp(itype(idih_constr(i)-1)),idih_constr(i)-1,
1009      &      phibound(1,i)*rad2deg,sdihed(1,i)*rad2deg,
1010      &      phibound(2,i)*rad2deg,sdihed(2,i)*rad2deg,
1011      &      (vpsipred(j,i),j=1,3)
1012          enddo
1013 #ifdef MPI
1014         endif
1015 #endif
1016 c Raw psipred input
1017       endif
1018 C first setting the theta boundaries to 0 to pi
1019 C this mean that there is no energy penalty for any angle occuring this can be applied
1020 C for generate random conformation but is not implemented in this way
1021 C      do i=1,nres
1022 C        thetabound(1,i)=0
1023 C        thetabound(2,i)=pi
1024 C      enddo
1025 C begin reading theta constrains this is quartic constrains allowing to
1026 C have smooth second derivative
1027       if (with_theta_constr) then
1028 C with_theta_constr is keyword allowing for occurance of theta constrains
1029       read (inp,*) ntheta_constr
1030 C ntheta_constr is the number of theta constrains
1031       if (ntheta_constr.gt.0) then
1032 C        read (inp,*) ftors
1033         read (inp,*) (itheta_constr(i),theta_constr0(i),
1034      &  theta_drange(i),for_thet_constr(i),
1035      &  i=1,ntheta_constr)
1036 C the above code reads from 1 to ntheta_constr
1037 C itheta_constr(i) residue i for which is theta_constr
1038 C theta_constr0 the global minimum value
1039 C theta_drange is range for which there is no energy penalty
1040 C for_thet_constr is the force constant for quartic energy penalty
1041 C E=k*x**4
1042         if(me.eq.king.or..not.out1file)then
1043          write (iout,*)
1044      &   'There are',ntheta_constr,' constraints on phi angles.'
1045          do i=1,ntheta_constr
1046           write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
1047      &    theta_drange(i),
1048      &    for_thet_constr(i)
1049          enddo
1050         endif
1051         do i=1,ntheta_constr
1052           theta_constr0(i)=deg2rad*theta_constr0(i)
1053           theta_drange(i)=deg2rad*theta_drange(i)
1054         enddo
1055 C        if(me.eq.king.or..not.out1file)
1056 C     &   write (iout,*) 'FTORS',ftors
1057 C        do i=1,ntheta_constr
1058 C          ii = itheta_constr(i)
1059 C          thetabound(1,ii) = phi0(i)-drange(i)
1060 C          thetabound(2,ii) = phi0(i)+drange(i)
1061 C        enddo
1062       endif ! ntheta_constr.gt.0
1063       endif! with_theta_constr
1064 C
1065 C      with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
1066 C      write (iout,*) "with_dihed_constr ",with_dihed_constr
1067 c--- Zscore rms -------
1068       if (nz_start.eq.0) nz_start=nnt
1069       if (nz_end.eq.0 .and. nsup.gt.0) then
1070         nz_end=nnt+nsup-1
1071       else if (nz_end.eq.0) then
1072         nz_end=nct
1073       endif
1074       if(me.eq.king.or..not.out1file)then
1075        write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1076        write (iout,*) 'IZ_SC=',iz_sc
1077       endif
1078 c----------------------
1079       call init_int_table
1080       if (refstr) then
1081         if (.not.pdbref) then
1082           call read_angles(inp,*38)
1083           call bond_regular
1084           goto 39
1085    38     write (iout,'(a)') 'Error reading reference structure.'
1086 #ifdef MPI
1087           call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1088           stop 'Error reading reference structure'
1089 #endif
1090    39     call chainbuild_extconf
1091           call setup_var
1092 czscore          call geom_to_var(nvar,coord_exp_zs(1,1))
1093           nstart_sup=nnt
1094           nstart_seq=nnt
1095           nsup=nct-nnt+1
1096           do i=1,2*nres
1097             do j=1,3
1098               cref(j,i)=c(j,i)
1099             enddo
1100           enddo
1101           call contact(.true.,ncont_ref,icont_ref,co)
1102         endif
1103       endif
1104 c        print *, "A TU"
1105 c      write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1106       call flush(iout)
1107       if (constr_dist.gt.0) call read_dist_constr
1108 c      write (iout,*) "After read_dist_constr nhpb",nhpb
1109       if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
1110       call hpb_partition
1111       call NMRpeak_partition
1112       if(me.eq.king.or..not.out1file)write (iout,*) 'Contact order:',co
1113       if (pdbref) then
1114         if(me.eq.king.or..not.out1file)
1115      &   write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1116         do i=1,ncont_ref
1117           do j=1,2
1118             icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1119           enddo
1120           if(me.eq.king.or..not.out1file)
1121      &     write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1122      &     icont_ref(1,i),' ',
1123      &     restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1124         enddo
1125       endif
1126       write (iout,*) "calling read_saxs_consrtr",nsaxs
1127       if (nsaxs.gt.0) call read_saxs_constr
1128
1129       if (constr_homology.gt.0) then
1130         call read_constr_homology
1131         if (indpdb.gt.0 .or. pdbref) then
1132           do i=1,2*nres
1133             do j=1,3
1134               c(j,i)=crefjlee(j,i)
1135               cref(j,i)=crefjlee(j,i)
1136             enddo
1137           enddo
1138         endif
1139 #ifdef DEBUG
1140         write (iout,*) "sc_loc_geom: Array C"
1141         do i=1,nres
1142           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),
1143      &      (c(j,i+nres),j=1,3)
1144         enddo
1145         write (iout,*) "Array Cref"
1146         do i=1,nres
1147           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i),j=1,3),
1148      &      (cref(j,i+nres),j=1,3)
1149         enddo
1150 #endif
1151        call int_from_cart1(.false.)
1152        call sc_loc_geom(.false.)
1153        do i=1,nres
1154          thetaref(i)=theta(i)
1155          phiref(i)=phi(i)
1156        enddo
1157        do i=1,nres-1
1158          do j=1,3
1159            dc(j,i)=c(j,i+1)-c(j,i)
1160            dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1161          enddo
1162        enddo
1163        do i=2,nres-1
1164          do j=1,3
1165            dc(j,i+nres)=c(j,i+nres)-c(j,i)
1166            dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1167          enddo
1168        enddo
1169       else
1170         homol_nset=0
1171         if (start_from_model) then
1172           nmodel_start=0
1173           do
1174             read(inp,'(a)',end=332,err=332) pdbfile
1175             if (me.eq.king .or. .not. out1file)
1176      &        write (iout,'(a,5x,a)') 'Opening PDB file',
1177      &        pdbfile(:ilen(pdbfile))
1178             open(ipdbin,file=pdbfile,status='old',err=336)
1179             goto 335
1180  336        write (iout,'(a,5x,a)') 'Error opening PDB file',
1181      &      pdbfile(:ilen(pdbfile))
1182             call flush(iout)
1183             stop
1184  335        continue
1185             unres_pdb=.false.
1186             nres_temp=nres
1187             call readpdb
1188             close(ipdbin)
1189             if (nres.ge.nres_temp) then
1190               nmodel_start=nmodel_start+1
1191               pdbfiles_chomo(nmodel_start)=pdbfile
1192               do i=1,2*nres
1193                 do j=1,3
1194                   chomo(j,i,nmodel_start)=c(j,i)
1195                 enddo
1196               enddo
1197             else
1198 c              itemp=nres
1199 c              nres=nres_temp
1200 c              call gen_rand_conf(itemp,*115)
1201 c              nmodel_start=nmodel_start+1
1202 c              do i=1,2*nres
1203 c                do j=1,3
1204 c                  chomo(j,i,nmodel_start)=c(j,i)
1205 c                enddo
1206 c              enddo
1207 c              goto 116
1208   115         if (me.eq.king .or. .not. out1file)
1209      &          write (iout,'(a,2i5,1x,a)')
1210      &           "Different number of residues",nres_temp,nres,
1211      &           " model skipped."
1212             endif
1213   116       continue
1214             nres=nres_temp
1215           enddo
1216   332     continue
1217           if (nmodel_start.eq.0) then
1218             if (me.eq.king .or. .not. out1file)
1219      &        write (iout,'(a)')
1220      &        "No valid starting model found START_FROM_MODELS is OFF"
1221               start_from_model=.false.
1222           endif
1223           write (iout,*) "nmodel_start",nmodel_start
1224         endif
1225       endif
1226
1227
1228 C      endif
1229       if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1230      &    .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1231      &    modecalc.ne.10) then
1232 C If input structure hasn't been supplied from the PDB file read or generate
1233 C initial geometry.
1234         if (iranconf.eq.0 .and. .not. extconf .and. .not.
1235      &     start_from_model) then
1236           if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1237      &     write (iout,'(a)') 'Initial geometry will be read in.'
1238           if (read_cart) then
1239             read(inp,'(8f10.5)',end=36,err=36)
1240      &       ((c(l,k),l=1,3),k=1,nres),
1241      &       ((c(l,k+nres),l=1,3),k=nnt,nct)
1242             if (nnt.gt.1) c(:,nres+1)=c(:,1)
1243             if (nct.lt.nres) c(:,2*nres)=c(:,nres)
1244 c            write (iout,*) "Exit READ_CART"
1245 c            write (iout,'(8f10.5)')
1246 c     &       ((c(l,k),l=1,3),k=1,nres),
1247 c     &       ((c(l,k+nres),l=1,3),k=nnt,nct)
1248             call cartprint
1249             do j=1,3
1250               dc(j,0)=c(j,1)
1251             enddo
1252             do i=1,nres-1
1253               do j=1,3
1254                 dc(j,i)=c(j,i+1)-c(j,i)
1255 c                dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1256               enddo
1257             enddo
1258             do i=nnt,nct
1259               if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1260                 do j=1,3
1261                   dc(j,i+nres)=c(j,i+nres)-c(j,i)
1262 c                  dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1263                 enddo
1264               else
1265                 do j=1,3
1266                   dc(j,i+nres)=0.0d0
1267 c                  dc_norm(j,i+nres)=0.0d0
1268                 enddo
1269               endif
1270             enddo
1271             call int_from_cart1(.true.)
1272             write (iout,*) "Finish INT_TO_CART"
1273 c      write (iout,*) "DC"
1274 c      do i=1,nres
1275 c        write (iout,'(i5,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1276 c     &    (dc(j,i+nres),j=1,3)
1277 c      enddo
1278 c            call cartprint
1279 c            call setup_var
1280 c            return
1281           else
1282             call read_angles(inp,*36)
1283             call bond_regular
1284             call chainbuild_extconf
1285           endif
1286           goto 37
1287    36     write (iout,'(a)') 'Error reading angle file.'
1288 #ifdef MPI
1289           call mpi_finalize( MPI_COMM_WORLD,IERR )
1290 #endif
1291           stop 'Error reading angle file.'
1292    37     continue
1293         else if (extconf) then
1294           if (me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1295      &      write (iout,'(a)') 'Extended chain initial geometry.'
1296           do i=3,nres
1297             theta(i)=90d0*deg2rad
1298           enddo
1299           do i=4,nres
1300             phi(i)=180d0*deg2rad
1301           enddo
1302           do i=2,nres-1
1303             alph(i)=110d0*deg2rad
1304           enddo
1305           do i=2,nres-1
1306             omeg(i)=-120d0*deg2rad
1307             if (itype(i).le.0) omeg(i)=-omeg(i)
1308           enddo
1309           call bond_regular
1310           call chainbuild_extconf
1311         else if (.not. start_from_model) then
1312           if(me.eq.king.or..not.out1file)
1313      &     write (iout,'(a)') 'Random-generated initial geometry.'
1314           call bond_regular
1315 #ifdef MPI
1316           if (me.eq.king  .or. fg_rank.eq.0 .and. (
1317      &           modecalc.eq.12 .or. modecalc.eq.14) ) then
1318 #endif
1319             do itrial=1,100
1320               itmp=1
1321               call gen_rand_conf(itmp,*30)
1322               goto 40
1323    30         write (iout,*) 'Failed to generate random conformation',
1324      &          ', itrial=',itrial
1325               write (*,*) 'Processor:',me,
1326      &          ' Failed to generate random conformation',
1327      &          ' itrial=',itrial
1328               call intout
1329
1330 #ifdef AIX
1331               call flush_(iout)
1332 #else
1333               call flush(iout)
1334 #endif
1335             enddo
1336             write (iout,'(a,i3,a)') 'Processor:',me,
1337      &        ' error in generating random conformation.'
1338             write (*,'(a,i3,a)') 'Processor:',me,
1339      &        ' error in generating random conformation.'
1340             call flush(iout)
1341 #ifdef MPI
1342             call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1343    40       continue
1344           endif
1345 #else
1346           write (*,'(a)')
1347      &      ' error in generating random conformation.'
1348           stop
1349    40     continue
1350 #endif
1351         endif
1352       elseif (modecalc.eq.4) then
1353         read (inp,'(a)') intinname
1354         open (intin,file=intinname,status='old',err=333)
1355         if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1356      &  write (iout,'(a)') 'intinname',intinname
1357         write (*,'(a)') 'Processor',myrank,' intinname',intinname
1358         goto 334
1359   333   write (iout,'(2a)') 'Error opening angle file ',intinname
1360 #ifdef MPI
1361         call MPI_Finalize(MPI_COMM_WORLD,IERR)
1362 #endif
1363         stop 'Error opening angle file.'
1364   334   continue
1365
1366       endif
1367 C Generate distance constraints, if the PDB structure is to be regularized.
1368       if (nthread.gt.0) then
1369         call read_threadbase
1370       endif
1371       call setup_var
1372       if (me.eq.king .or. .not. out1file)
1373      & call intout
1374       if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1375         write (iout,'(/a,i3,a)')
1376      &  'The chain contains',ns,' disulfide-bridging cysteines.'
1377         write (iout,'(20i4)') (iss(i),i=1,ns)
1378        if (dyn_ss) then
1379           write(iout,*)"Running with dynamic disulfide-bond formation"
1380        else
1381         write (iout,'(/a/)') 'Pre-formed links are:'
1382         do i=1,nss
1383           i1=ihpb(i)-nres
1384           i2=jhpb(i)-nres
1385           it1=itype(i1)
1386           it2=itype(i2)
1387           write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1388      &    restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1389      &    ebr,forcon(i)
1390         enddo
1391         write (iout,'(a)')
1392        endif
1393       endif
1394       if (ns.gt.0.and.dyn_ss) then
1395           do i=nss+1,nhpb
1396             ihpb(i-nss)=ihpb(i)
1397             jhpb(i-nss)=jhpb(i)
1398             forcon(i-nss)=forcon(i)
1399             dhpb(i-nss)=dhpb(i)
1400           enddo
1401           nhpb=nhpb-nss
1402           nss=0
1403           call hpb_partition
1404           do i=1,ns
1405             dyn_ss_mask(iss(i))=.true.
1406           enddo
1407       endif
1408 c --- test
1409 c      call cartprint
1410 c      write (iout,*) "DC"
1411 c      do i=1,nres
1412 c        write (iout,'(i5,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1413 c     &    (dc(j,i+nres),j=1,3)
1414 c      enddo
1415       if (i2ndstr.gt.0) call secstrp2dihc
1416 c      call geom_to_var(nvar,x)
1417 c      call etotal(energia(0))
1418 c      call enerprint(energia(0))
1419 c      call briefout(0,etot)
1420 c      stop
1421 cd    write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1422 cd    write (iout,'(a)') 'Variable list:'
1423 cd    write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1424 #ifdef MPI
1425       if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1426      &  write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1427      &  'Processor',myrank,': end reading molecular data.'
1428 #endif
1429 c      print *,"A TU?"
1430       return
1431       end
1432 c--------------------------------------------------------------------------
1433       logical function seq_comp(itypea,itypeb,length)
1434       implicit none
1435       integer length,itypea(length),itypeb(length)
1436       integer i
1437       do i=1,length
1438         if (itypea(i).ne.itypeb(i)) then
1439           seq_comp=.false.
1440           return
1441         endif
1442       enddo
1443       seq_comp=.true.
1444       return
1445       end
1446 c-----------------------------------------------------------------------------
1447       subroutine read_bridge
1448 C Read information about disulfide bridges.
1449       implicit none
1450       include 'DIMENSIONS'
1451 #ifdef MPI
1452       include 'mpif.h'
1453       integer ierror
1454 #endif
1455       include 'COMMON.IOUNITS'
1456       include 'COMMON.GEO'
1457       include 'COMMON.VAR'
1458       include 'COMMON.INTERACT'
1459       include 'COMMON.LOCAL'
1460       include 'COMMON.NAMES'
1461       include 'COMMON.CHAIN'
1462       include 'COMMON.FFIELD'
1463       include 'COMMON.SBRIDGE'
1464       include 'COMMON.HEADER'
1465       include 'COMMON.CONTROL'
1466       include 'COMMON.DBASE'
1467       include 'COMMON.THREAD'
1468       include 'COMMON.TIME1'
1469       include 'COMMON.SETUP'
1470       integer i,j
1471 C Read bridging residues.
1472       read (inp,*) ns,(iss(i),i=1,ns)
1473 c 5/24/2020 Adam: Added a table to translate residue numbers to cysteine
1474 c     numbers
1475       icys=0
1476       do i=1,ns
1477         icys(iss(i))=i
1478       enddo
1479 c      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 (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 c          dhpb(i)=dbr
1533 c          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 none
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       integer i,j,k,l,kanal
1555 c Read coordinates from input
1556 c
1557       read(kanal,'(8f10.5)',end=10,err=10)
1558      &  ((c(l,k),l=1,3),k=1,nres),
1559      &  ((c(l,k+nres),l=1,3),k=nnt,nct)
1560       do j=1,3
1561         c(j,nres+1)=c(j,1)
1562         c(j,2*nres)=c(j,nres)
1563       enddo
1564       call int_from_cart1(.false.)
1565       do i=1,nres-1
1566         do j=1,3
1567           dc(j,i)=c(j,i+1)-c(j,i)
1568           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1569         enddo
1570       enddo
1571       do i=nnt,nct
1572         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1573           do j=1,3
1574             dc(j,i+nres)=c(j,i+nres)-c(j,i)
1575             dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1576           enddo
1577         endif
1578       enddo
1579
1580       return
1581    10 return1
1582       end
1583 c----------------------------------------------------------------------------
1584       subroutine read_threadbase
1585       implicit none
1586       include 'DIMENSIONS'
1587       include 'COMMON.IOUNITS'
1588       include 'COMMON.GEO'
1589       include 'COMMON.VAR'
1590       include 'COMMON.INTERACT'
1591       include 'COMMON.LOCAL'
1592       include 'COMMON.NAMES'
1593       include 'COMMON.CHAIN'
1594       include 'COMMON.FFIELD'
1595       include 'COMMON.SBRIDGE'
1596       include 'COMMON.HEADER'
1597       include 'COMMON.CONTROL'
1598       include 'COMMON.DBASE'
1599       include 'COMMON.THREAD'
1600       include 'COMMON.TIME1'
1601       integer i,j,k
1602       double precision dist
1603 C Read pattern database for threading.
1604       read (icbase,*) nseq
1605       do i=1,nseq
1606         read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1607      &   nres_base(2,i),nres_base(3,i)
1608         read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1609      &   nres_base(1,i))
1610 c       write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1611 c    &   nres_base(2,i),nres_base(3,i)
1612 c       write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1613 c    &   nres_base(1,i))
1614       enddo
1615       close (icbase)
1616       if (weidis.eq.0.0D0) weidis=0.1D0
1617       do i=nnt,nct
1618         do j=i+2,nct
1619           nhpb=nhpb+1
1620           ihpb(nhpb)=i
1621           jhpb(nhpb)=j
1622           forcon(nhpb)=weidis
1623         enddo
1624       enddo
1625       read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1626       write (iout,'(a,i5)') 'nexcl: ',nexcl
1627       write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1628       return
1629       end
1630 c------------------------------------------------------------------------------
1631       subroutine setup_var
1632       implicit none
1633       include 'DIMENSIONS'
1634       include 'COMMON.IOUNITS'
1635       include 'COMMON.GEO'
1636       include 'COMMON.VAR'
1637       include 'COMMON.INTERACT'
1638       include 'COMMON.LOCAL'
1639       include 'COMMON.NAMES'
1640       include 'COMMON.CHAIN'
1641       include 'COMMON.FFIELD'
1642       include 'COMMON.SBRIDGE'
1643       include 'COMMON.HEADER'
1644       include 'COMMON.CONTROL'
1645       include 'COMMON.DBASE'
1646       include 'COMMON.THREAD'
1647       include 'COMMON.TIME1'
1648       integer i
1649 C Set up variable list.
1650       ntheta=nres-2
1651       nphi=nres-3
1652       nvar=ntheta+nphi
1653       nside=0
1654       do i=2,nres-1
1655         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1656           nside=nside+1
1657           ialph(i,1)=nvar+nside
1658           ialph(nside,2)=i
1659         endif
1660       enddo
1661       if (indphi.gt.0) then
1662         nvar=nphi
1663       else if (indback.gt.0) then
1664         nvar=nphi+ntheta
1665       else
1666         nvar=nvar+2*nside
1667       endif
1668       return
1669       end
1670 c----------------------------------------------------------------------------
1671       subroutine gen_dist_constr
1672 C Generate CA distance constraints.
1673       implicit none
1674       include 'DIMENSIONS'
1675       include 'COMMON.IOUNITS'
1676       include 'COMMON.GEO'
1677       include 'COMMON.VAR'
1678       include 'COMMON.INTERACT'
1679       include 'COMMON.LOCAL'
1680       include 'COMMON.NAMES'
1681       include 'COMMON.CHAIN'
1682       include 'COMMON.FFIELD'
1683       include 'COMMON.SBRIDGE'
1684       include 'COMMON.HEADER'
1685       include 'COMMON.CONTROL'
1686       include 'COMMON.DBASE'
1687       include 'COMMON.THREAD'
1688       include 'COMMON.SPLITELE'
1689       include 'COMMON.TIME1'
1690       integer i,j,itype_pdb(maxres)
1691       common /pizda/ itype_pdb
1692       double precision dd
1693       double precision dist
1694       character*2 iden
1695 cd      print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1696 cd      write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1697 cd     & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1698 cd     & ' nsup',nsup
1699       do i=nstart_sup,nstart_sup+nsup-1
1700 cd      write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1701 cd     &    ' seq_pdb', restyp(itype_pdb(i))
1702         do j=i+2,nstart_sup+nsup-1
1703 c 5/24/2020 Adam: Cutoff included to reduce array size
1704           dd = dist(i,j)
1705           if (dd.gt.r_cut_int) cycle
1706           nhpb=nhpb+1
1707           ihpb(nhpb)=i+nstart_seq-nstart_sup
1708           jhpb(nhpb)=j+nstart_seq-nstart_sup
1709           forcon(nhpb)=weidis
1710           dhpb(nhpb)=dd
1711         enddo
1712       enddo
1713 cd      write (iout,'(a)') 'Distance constraints:'
1714 cd      do i=nss+1,nhpb
1715 cd        ii=ihpb(i)
1716 cd        jj=jhpb(i)
1717 cd        iden='CA'
1718 cd        if (ii.gt.nres) then
1719 cd          iden='SC'
1720 cd          ii=ii-nres
1721 cd          jj=jj-nres
1722 cd        endif
1723 cd        write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1724 cd     &  restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1725 cd     &  dhpb(i),forcon(i)
1726 cd      enddo
1727       return
1728       end
1729 c----------------------------------------------------------------------------
1730       subroutine map_read
1731       implicit none
1732       include 'DIMENSIONS'
1733       include 'COMMON.MAP'
1734       include 'COMMON.IOUNITS'
1735       integer imap
1736       character*3 angid(4) /'THE','PHI','ALP','OME'/
1737       character*80 mapcard,ucase
1738       do imap=1,nmap
1739         read (inp,'(a)') mapcard
1740         mapcard=ucase(mapcard)
1741         if (index(mapcard,'PHI').gt.0) then
1742           kang(imap)=1
1743         else if (index(mapcard,'THE').gt.0) then
1744           kang(imap)=2
1745         else if (index(mapcard,'ALP').gt.0) then
1746           kang(imap)=3
1747         else if (index(mapcard,'OME').gt.0) then
1748           kang(imap)=4
1749         else
1750           write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1751           stop 'Error - illegal variable spec in MAP card.'
1752         endif
1753         call readi (mapcard,'RES1',res1(imap),0)
1754         call readi (mapcard,'RES2',res2(imap),0)
1755         if (res1(imap).eq.0) then
1756           res1(imap)=res2(imap)
1757         else if (res2(imap).eq.0) then
1758           res2(imap)=res1(imap)
1759         endif
1760         if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1761           write (iout,'(a)')
1762      &    'Error - illegal definition of variable group in MAP.'
1763           stop 'Error - illegal definition of variable group in MAP.'
1764         endif
1765         call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1766         call reada(mapcard,'TO',ang_to(imap),0.0D0)
1767         call readi(mapcard,'NSTEP',nstep(imap),0)
1768         if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1769           write (iout,'(a)')
1770      &     'Illegal boundary and/or step size specification in MAP.'
1771           stop 'Illegal boundary and/or step size specification in MAP.'
1772         endif
1773       enddo ! imap
1774       return
1775       end
1776 c----------------------------------------------------------------------------
1777       subroutine csaread
1778       implicit none
1779       include 'DIMENSIONS'
1780       include 'COMMON.IOUNITS'
1781       include 'COMMON.GEO'
1782       include 'COMMON.CSA'
1783       include 'COMMON.BANK'
1784       include 'COMMON.CONTROL'
1785       character*80 ucase
1786       character*620 mcmcard
1787       call card_concat(mcmcard)
1788
1789       call readi(mcmcard,'NCONF',nconf,50)
1790       call readi(mcmcard,'NADD',nadd,0)
1791       call readi(mcmcard,'JSTART',jstart,1)
1792       call readi(mcmcard,'JEND',jend,1)
1793       call readi(mcmcard,'NSTMAX',nstmax,500000)
1794       call readi(mcmcard,'N0',n0,1)
1795       call readi(mcmcard,'N1',n1,6)
1796       call readi(mcmcard,'N2',n2,4)
1797       call readi(mcmcard,'N3',n3,0)
1798       call readi(mcmcard,'N4',n4,0)
1799       call readi(mcmcard,'N5',n5,0)
1800       call readi(mcmcard,'N6',n6,10)
1801       call readi(mcmcard,'N7',n7,0)
1802       call readi(mcmcard,'N8',n8,0)
1803       call readi(mcmcard,'N9',n9,0)
1804       call readi(mcmcard,'N14',n14,0)
1805       call readi(mcmcard,'N15',n15,0)
1806       call readi(mcmcard,'N16',n16,0)
1807       call readi(mcmcard,'N17',n17,0)
1808       call readi(mcmcard,'N18',n18,0)
1809
1810       vdisulf=(index(mcmcard,'DYNSS').gt.0)
1811
1812       call readi(mcmcard,'NDIFF',ndiff,2)
1813       call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1814       call readi(mcmcard,'IS1',is1,1)
1815       call readi(mcmcard,'IS2',is2,8)
1816       call readi(mcmcard,'NRAN0',nran0,4)
1817       call readi(mcmcard,'NRAN1',nran1,2)
1818       call readi(mcmcard,'IRR',irr,1)
1819       call readi(mcmcard,'NSEED',nseed,20)
1820       call readi(mcmcard,'NTOTAL',ntotal,10000)
1821       call reada(mcmcard,'CUT1',cut1,2.0d0)
1822       call reada(mcmcard,'CUT2',cut2,5.0d0)
1823       call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1824       call readi(mcmcard,'ICMAX',icmax,3)
1825       call readi(mcmcard,'IRESTART',irestart,0)
1826 c!bankt      call readi(mcmcard,'NBANKTM',ntbankm,0)
1827       ntbankm=0
1828 c!bankt
1829       call reada(mcmcard,'DELE',dele,20.0d0)
1830       call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1831       call readi(mcmcard,'IREF',iref,0)
1832       call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1833       call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1834       call readi(mcmcard,'NCONF_IN',nconf_in,0)
1835       call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1836       write (iout,*) "NCONF_IN",nconf_in
1837       return
1838       end
1839 c----------------------------------------------------------------------------
1840       subroutine mcmread
1841       implicit none
1842       include 'DIMENSIONS'
1843       include 'COMMON.MCM'
1844       include 'COMMON.MCE'
1845       include 'COMMON.IOUNITS'
1846       character*80 ucase
1847       character*320 mcmcard
1848       integer i
1849       call card_concat(mcmcard)
1850       call readi(mcmcard,'MAXACC',maxacc,100)
1851       call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1852       call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1853       call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1854       call readi(mcmcard,'MAXREPM',maxrepm,200)
1855       call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1856       call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1857       call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1858       call reada(mcmcard,'E_UP',e_up,5.0D0)
1859       call reada(mcmcard,'DELTE',delte,0.1D0)
1860       call readi(mcmcard,'NSWEEP',nsweep,5)
1861       call readi(mcmcard,'NSTEPH',nsteph,0)
1862       call readi(mcmcard,'NSTEPC',nstepc,0)
1863       call reada(mcmcard,'TMIN',tmin,298.0D0)
1864       call reada(mcmcard,'TMAX',tmax,298.0D0)
1865       call readi(mcmcard,'NWINDOW',nwindow,0)
1866       call readi(mcmcard,'PRINT_MC',print_mc,0)
1867       print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1868       print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1869       ent_read=(index(mcmcard,'ENT_READ').gt.0)
1870       call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1871       call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1872       call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1873       call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1874       call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1875       if (nwindow.gt.0) then
1876         read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1877         do i=1,nwindow
1878           winlen(i)=winend(i)-winstart(i)+1
1879         enddo
1880       endif
1881       if (tmax.lt.tmin) tmax=tmin
1882       if (tmax.eq.tmin) then
1883         nstepc=0
1884         nsteph=0
1885       endif
1886       if (nstepc.gt.0 .and. nsteph.gt.0) then
1887         tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1888         tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1889       endif
1890 C Probabilities of different move types
1891       sumpro_type(0)=0.0D0
1892       call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1893       call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1894       sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1895       call reada(mcmcard,'THETA'     ,sumpro_type(3),0.0d0)
1896       sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1897       call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1898       sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1899       do i=1,MaxMoveType
1900         print *,'i',i,' sumprotype',sumpro_type(i)
1901         sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1902         print *,'i',i,' sumprotype',sumpro_type(i)
1903       enddo
1904       return
1905       end
1906 c----------------------------------------------------------------------------
1907       subroutine read_minim
1908       implicit none
1909       include 'DIMENSIONS'
1910       include 'COMMON.MINIM'
1911       include 'COMMON.IOUNITS'
1912       character*80 ucase
1913       character*320 minimcard
1914       call card_concat(minimcard)
1915       call readi(minimcard,'MAXMIN',maxmin,2000)
1916       call readi(minimcard,'MAXFUN',maxfun,5000)
1917       call readi(minimcard,'MINMIN',minmin,maxmin)
1918       call readi(minimcard,'MINFUN',minfun,maxmin)
1919       call reada(minimcard,'TOLF',tolf,1.0D-2)
1920       call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1921       print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1922       print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1923       print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1924       write (iout,'(/80(1h*)/20x,a/80(1h*))')
1925      &         'Options in energy minimization:'
1926       write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1927      & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1928      & 'MinMin:',MinMin,' MinFun:',MinFun,
1929      & ' TolF:',TolF,' RTolF:',RTolF
1930       return
1931       end
1932 c----------------------------------------------------------------------------
1933       subroutine read_angles(kanal,*)
1934       implicit none
1935       include 'DIMENSIONS'
1936       include 'COMMON.GEO'
1937       include 'COMMON.VAR'
1938       include 'COMMON.CHAIN'
1939       include 'COMMON.IOUNITS'
1940       include 'COMMON.CONTROL'
1941       integer i,kanal
1942 c Read angles from input
1943 c
1944        read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1945        read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1946        read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1947        read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1948
1949        do i=1,nres
1950 c 9/7/01 avoid 180 deg valence angle
1951         if (theta(i).gt.179.99d0) theta(i)=179.99d0
1952 c
1953         theta(i)=deg2rad*theta(i)
1954         phi(i)=deg2rad*phi(i)
1955         alph(i)=deg2rad*alph(i)
1956         omeg(i)=deg2rad*omeg(i)
1957        enddo
1958       return
1959    10 return1
1960       end
1961 c----------------------------------------------------------------------------
1962       subroutine reada(rekord,lancuch,wartosc,default)
1963       implicit none
1964       character*(*) rekord,lancuch
1965       double precision wartosc,default
1966       integer ilen,iread
1967       external ilen
1968       iread=index(rekord,lancuch)
1969       if (iread.eq.0) then
1970         wartosc=default
1971         return
1972       endif
1973       iread=iread+ilen(lancuch)+1
1974       read (rekord(iread:),*,err=10,end=10) wartosc
1975       return
1976   10  wartosc=default
1977       return
1978       end
1979 c----------------------------------------------------------------------------
1980       subroutine readi(rekord,lancuch,wartosc,default)
1981       implicit none
1982       character*(*) rekord,lancuch
1983       integer wartosc,default
1984       integer ilen,iread
1985       external ilen
1986       iread=index(rekord,lancuch)
1987       if (iread.eq.0) then
1988         wartosc=default
1989         return
1990       endif
1991       iread=iread+ilen(lancuch)+1
1992       read (rekord(iread:),*,err=10,end=10) wartosc
1993       return
1994   10  wartosc=default
1995       return
1996       end
1997 c----------------------------------------------------------------------------
1998       subroutine multreadi(rekord,lancuch,tablica,dim,default)
1999       implicit none
2000       integer dim,i
2001       integer tablica(dim),default
2002       character*(*) rekord,lancuch
2003       character*80 aux
2004       integer ilen,iread
2005       external ilen
2006       do i=1,dim
2007         tablica(i)=default
2008       enddo
2009       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2010       if (iread.eq.0) return
2011       iread=iread+ilen(lancuch)+1
2012       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2013    10 return
2014       end
2015 c----------------------------------------------------------------------------
2016       subroutine multreada(rekord,lancuch,tablica,dim,default)
2017       implicit none
2018       integer dim,i
2019       double precision tablica(dim),default
2020       character*(*) rekord,lancuch
2021       character*80 aux
2022       integer ilen,iread
2023       external ilen
2024       do i=1,dim
2025         tablica(i)=default
2026       enddo
2027       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2028       if (iread.eq.0) return
2029       iread=iread+ilen(lancuch)+1
2030       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2031    10 return
2032       end
2033 c----------------------------------------------------------------------------
2034       subroutine openunits
2035       implicit none
2036       include 'DIMENSIONS'
2037 #ifdef MPI
2038       include 'mpif.h'
2039       integer ierror
2040       character*16 form,nodename
2041       integer nodelen
2042 #endif
2043       include 'COMMON.SETUP'
2044       include 'COMMON.IOUNITS'
2045       include 'COMMON.MD'
2046       include 'COMMON.CONTROL'
2047       integer lenpre,lenpot,ilen,lentmp,npos
2048       external ilen
2049       character*3 out1file_text,ucase
2050       character*3 ll
2051       external ucase
2052       tmpdir=""
2053       out1file_text="YES"
2054 c      print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
2055       call getenv_loc("PREFIX",prefix)
2056       pref_orig = prefix
2057       call getenv_loc("POT",pot)
2058       call getenv_loc("DIRTMP",tmpdir)
2059       call getenv_loc("CURDIR",curdir)
2060       call getenv_loc("OUT1FILE",out1file_text)
2061 c      print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
2062       out1file_text=ucase(out1file_text)
2063       if (out1file_text(1:1).eq."Y") then
2064         out1file=.true.
2065       else
2066         out1file=fg_rank.gt.0
2067       endif
2068       lenpre=ilen(prefix)
2069       lenpot=ilen(pot)
2070       lentmp=ilen(tmpdir)
2071       if (lentmp.gt.0) then
2072           write (*,'(80(1h!))')
2073           write (*,'(a,19x,a,19x,a)') "!","  A T T E N T I O N  ","!"
2074           write (*,'(80(1h!))')
2075           write (*,*)"All output files will be on node /tmp directory."
2076 #ifdef MPI
2077         call  MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
2078         if (me.eq.king) then
2079           write (*,*) "The master node is ",nodename
2080         else if (fg_rank.eq.0) then
2081           write (*,*) "I am the CG slave node ",nodename
2082         else
2083           write (*,*) "I am the FG slave node ",nodename
2084         endif
2085 #endif
2086         PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
2087         lenpre = lentmp+lenpre+1
2088       endif
2089       entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
2090 C Get the names and open the input files
2091 #if defined(WINIFL) || defined(WINPGI)
2092       open(1,file=pref_orig(:ilen(pref_orig))//
2093      &  '.inp',status='old',readonly,shared)
2094        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2095 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2096 C Get parameter filenames and open the parameter files.
2097       call getenv_loc('BONDPAR',bondname)
2098       open (ibond,file=bondname,status='old',readonly,shared)
2099       call getenv_loc('THETPAR',thetname)
2100       open (ithep,file=thetname,status='old',readonly,shared)
2101       call getenv_loc('ROTPAR',rotname)
2102       open (irotam,file=rotname,status='old',readonly,shared)
2103       call getenv_loc('TORPAR',torname)
2104       open (itorp,file=torname,status='old',readonly,shared)
2105       call getenv_loc('TORDPAR',tordname)
2106       open (itordp,file=tordname,status='old',readonly,shared)
2107       call getenv_loc('FOURIER',fouriername)
2108       open (ifourier,file=fouriername,status='old',readonly,shared)
2109       call getenv_loc('ELEPAR',elename)
2110       open (ielep,file=elename,status='old',readonly,shared)
2111       call getenv_loc('SIDEPAR',sidename)
2112       open (isidep,file=sidename,status='old',readonly,shared)
2113       call getenv_loc('LIPTRANPAR',liptranname)
2114       open (iliptranpar,file=liptranname,status='old',readonly,shared)
2115 #elif (defined CRAY) || (defined AIX)
2116       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2117      &  action='read')
2118 c      print *,"Processor",myrank," opened file 1"
2119       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2120 c      print *,"Processor",myrank," opened file 9"
2121 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2122 C Get parameter filenames and open the parameter files.
2123       call getenv_loc('BONDPAR',bondname)
2124       open (ibond,file=bondname,status='old',action='read')
2125 c      print *,"Processor",myrank," opened file IBOND"
2126       call getenv_loc('THETPAR',thetname)
2127       open (ithep,file=thetname,status='old',action='read')
2128 #ifndef CRYST_THETA
2129       call getenv_loc('THETPARPDB',thetname_pdb)
2130       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2131 #endif
2132 c      print *,"Processor",myrank," opened file ITHEP"
2133       call getenv_loc('ROTPAR',rotname)
2134       open (irotam,file=rotname,status='old',action='read')
2135 #ifndef CRYST_SC
2136       call getenv_loc('ROTPARPDB',rotname_pdb)
2137       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2138 #endif
2139 c      print *,"Processor",myrank," opened file IROTAM"
2140       call getenv_loc('TORPAR',torname)
2141       open (itorp,file=torname,status='old',action='read')
2142 c      print *,"Processor",myrank," opened file ITORP"
2143       call getenv_loc('TORDPAR',tordname)
2144       open (itordp,file=tordname,status='old',action='read')
2145 c      print *,"Processor",myrank," opened file ITORDP"
2146       call getenv_loc('SCCORPAR',sccorname)
2147       open (isccor,file=sccorname,status='old',action='read')
2148 c      print *,"Processor",myrank," opened file ISCCOR"
2149       call getenv_loc('FOURIER',fouriername)
2150       open (ifourier,file=fouriername,status='old',action='read')
2151 c      print *,"Processor",myrank," opened file IFOURIER"
2152       call getenv_loc('ELEPAR',elename)
2153       open (ielep,file=elename,status='old',action='read')
2154 c      print *,"Processor",myrank," opened file IELEP"
2155       call getenv_loc('SIDEPAR',sidename)
2156       open (isidep,file=sidename,status='old',action='read')
2157       call getenv_loc('LIPTRANPAR',liptranname)
2158       open (iliptranpar,file=liptranname,status='old',action='read')
2159 c      print *,"Processor",myrank," opened file ISIDEP"
2160 c      print *,"Processor",myrank," opened parameter files"
2161 #elif (defined G77)
2162       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2163       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2164 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2165 C Get parameter filenames and open the parameter files.
2166       call getenv_loc('BONDPAR',bondname)
2167       open (ibond,file=bondname,status='old')
2168       call getenv_loc('THETPAR',thetname)
2169       open (ithep,file=thetname,status='old')
2170 #ifndef CRYST_THETA
2171       call getenv_loc('THETPARPDB',thetname_pdb)
2172       open (ithep_pdb,file=thetname_pdb,status='old')
2173 #endif
2174       call getenv_loc('ROTPAR',rotname)
2175       open (irotam,file=rotname,status='old')
2176 #ifndef CRYST_SC
2177       call getenv_loc('ROTPARPDB',rotname_pdb)
2178       open (irotam_pdb,file=rotname_pdb,status='old')
2179 #endif
2180       call getenv_loc('TORPAR',torname)
2181       open (itorp,file=torname,status='old')
2182       call getenv_loc('TORDPAR',tordname)
2183       open (itordp,file=tordname,status='old')
2184       call getenv_loc('SCCORPAR',sccorname)
2185       open (isccor,file=sccorname,status='old')
2186       call getenv_loc('FOURIER',fouriername)
2187       open (ifourier,file=fouriername,status='old')
2188       call getenv_loc('ELEPAR',elename)
2189       open (ielep,file=elename,status='old')
2190       call getenv_loc('SIDEPAR',sidename)
2191       open (isidep,file=sidename,status='old')
2192       call getenv_loc('LIPTRANPAR',liptranname)
2193       open (iliptranpar,file=liptranname,status='old')
2194 #else
2195       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2196      &  readonly)
2197        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2198 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2199 C Get parameter filenames and open the parameter files.
2200       call getenv_loc('BONDPAR',bondname)
2201       open (ibond,file=bondname,status='old',readonly)
2202       call getenv_loc('THETPAR',thetname)
2203       open (ithep,file=thetname,status='old',readonly)
2204       call getenv_loc('ROTPAR',rotname)
2205       open (irotam,file=rotname,status='old',readonly)
2206       call getenv_loc('TORPAR',torname)
2207       open (itorp,file=torname,status='old',readonly)
2208       call getenv_loc('TORDPAR',tordname)
2209       open (itordp,file=tordname,status='old',readonly)
2210       call getenv_loc('SCCORPAR',sccorname)
2211       open (isccor,file=sccorname,status='old',readonly)
2212 #ifndef CRYST_THETA
2213       call getenv_loc('THETPARPDB',thetname_pdb)
2214       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2215 #endif
2216       call getenv_loc('FOURIER',fouriername)
2217       open (ifourier,file=fouriername,status='old',readonly)
2218       call getenv_loc('ELEPAR',elename)
2219       open (ielep,file=elename,status='old',readonly)
2220       call getenv_loc('SIDEPAR',sidename)
2221       open (isidep,file=sidename,status='old',readonly)
2222       call getenv_loc('LIPTRANPAR',liptranname)
2223       open (iliptranpar,file=liptranname,status='old',action='read')
2224 #ifndef CRYST_SC
2225       call getenv_loc('ROTPARPDB',rotname_pdb)
2226       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2227 #endif
2228 #endif
2229 #ifdef TUBE
2230       call getenv_loc('TUBEPAR',tubename)
2231 #if defined(WINIFL) || defined(WINPGI)
2232       open (itube,file=tubename,status='old',readonly,shared)
2233 #elif (defined CRAY)  || (defined AIX)
2234       open (itube,file=tubename,status='old',action='read')
2235 #elif (defined G77)
2236       open (itube,file=tubename,status='old')
2237 #else
2238       open (itube,file=tubename,status='old',readonly)
2239 #endif
2240 #endif
2241 #ifndef OLDSCP
2242 C
2243 C 8/9/01 In the newest version SCp interaction constants are read from a file
2244 C Use -DOLDSCP to use hard-coded constants instead.
2245 C
2246       call getenv_loc('SCPPAR',scpname)
2247 #if defined(WINIFL) || defined(WINPGI)
2248       open (iscpp,file=scpname,status='old',readonly,shared)
2249 #elif (defined CRAY)  || (defined AIX)
2250       open (iscpp,file=scpname,status='old',action='read')
2251 #elif (defined G77)
2252       open (iscpp,file=scpname,status='old')
2253 #else
2254       open (iscpp,file=scpname,status='old',readonly)
2255 #endif
2256 #endif
2257       call getenv_loc('PATTERN',patname)
2258 #if defined(WINIFL) || defined(WINPGI)
2259       open (icbase,file=patname,status='old',readonly,shared)
2260 #elif (defined CRAY)  || (defined AIX)
2261       open (icbase,file=patname,status='old',action='read')
2262 #elif (defined G77)
2263       open (icbase,file=patname,status='old')
2264 #else
2265       open (icbase,file=patname,status='old',readonly)
2266 #endif
2267 #ifdef MPI
2268 C Open output file only for CG processes
2269 c      print *,"Processor",myrank," fg_rank",fg_rank
2270       if (fg_rank.eq.0) then
2271
2272       if (nodes.eq.1) then
2273         npos=3
2274       else
2275         npos = dlog10(dfloat(nodes-1))+1
2276       endif
2277       if (npos.lt.3) npos=3
2278       write (liczba,'(i1)') npos
2279       form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2280      &  //')'
2281       write (liczba,form) me
2282       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2283      &  liczba(:ilen(liczba))
2284       intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2285      &  //'.int'
2286       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2287      &  //'.pdb'
2288       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2289      &  liczba(:ilen(liczba))//'.mol2'
2290       statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2291      &  liczba(:ilen(liczba))//'.stat'
2292       if (lentmp.gt.0)
2293      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2294      &      //liczba(:ilen(liczba))//'.stat')
2295       rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2296      &  //'.rst'
2297       if(usampl) then
2298           qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2299      & liczba(:ilen(liczba))//'.const'
2300       endif
2301
2302       endif
2303 #else
2304       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2305       intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2306       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2307       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2308       statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2309       if (lentmp.gt.0)
2310      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2311      &    //'.stat')
2312       rest2name=prefix(:ilen(prefix))//'.rst'
2313       if(usampl) then
2314          qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2315       endif
2316 #endif
2317 #if defined(AIX) || defined(PGI) || defined(CRAY)
2318       if (me.eq.king .or. .not. out1file) then
2319          open(iout,file=outname,status='unknown')
2320       else
2321          open(iout,file="/dev/null",status="unknown")
2322       endif
2323 c#define DEBUG
2324 #ifdef DEBUG
2325       if (fg_rank.gt.0) then
2326         write (liczba,'(i3.3)') myrank/nfgtasks
2327         write (ll,'(bz,i3.3)') fg_rank
2328         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2329      &   status='unknown')
2330       endif
2331 #endif
2332 c#undef DEBUG
2333       if(me.eq.king) then
2334        open(igeom,file=intname,status='unknown',position='append')
2335        open(ipdb,file=pdbname,status='unknown')
2336        open(imol2,file=mol2name,status='unknown')
2337        open(istat,file=statname,status='unknown',position='append')
2338       else
2339 c1out       open(iout,file=outname,status='unknown')
2340       endif
2341 #else
2342       if (me.eq.king .or. .not.out1file)
2343      &    open(iout,file=outname,status='unknown')
2344 #ifdef DEBUG
2345       if (fg_rank.gt.0) then
2346         write (liczba,'(i3.3)') myrank/nfgtasks
2347         write (ll,'(bz,i3.3)') fg_rank
2348         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2349      &   status='unknown')
2350       endif
2351 #endif
2352       if(me.eq.king) then
2353        open(igeom,file=intname,status='unknown',access='append')
2354        open(ipdb,file=pdbname,status='unknown')
2355        open(imol2,file=mol2name,status='unknown')
2356        open(istat,file=statname,status='unknown',access='append')
2357       else
2358 c1out       open(iout,file=outname,status='unknown')
2359       endif
2360 #endif
2361       csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2362       csa_seed=prefix(:lenpre)//'.CSA.seed'
2363       csa_history=prefix(:lenpre)//'.CSA.history'
2364       csa_bank=prefix(:lenpre)//'.CSA.bank'
2365       csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2366       csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2367       csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2368 c!bankt      csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2369       csa_int=prefix(:lenpre)//'.int'
2370       csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2371       csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2372       csa_in=prefix(:lenpre)//'.CSA.in'
2373 c      print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2374 C Write file names
2375       if (me.eq.king)then
2376       write (iout,'(80(1h-))')
2377       write (iout,'(30x,a)') "FILE ASSIGNMENT"
2378       write (iout,'(80(1h-))')
2379       write (iout,*) "Input file                      : ",
2380      &  pref_orig(:ilen(pref_orig))//'.inp'
2381       write (iout,*) "Output file                     : ",
2382      &  outname(:ilen(outname))
2383       write (iout,*)
2384       write (iout,*) "Sidechain potential file        : ",
2385      &  sidename(:ilen(sidename))
2386 #ifndef OLDSCP
2387       write (iout,*) "SCp potential file              : ",
2388      &  scpname(:ilen(scpname))
2389 #endif
2390       write (iout,*) "Electrostatic potential file    : ",
2391      &  elename(:ilen(elename))
2392       write (iout,*) "Cumulant coefficient file       : ",
2393      &  fouriername(:ilen(fouriername))
2394       write (iout,*) "Torsional parameter file        : ",
2395      &  torname(:ilen(torname))
2396       write (iout,*) "Double torsional parameter file : ",
2397      &  tordname(:ilen(tordname))
2398       write (iout,*) "SCCOR parameter file : ",
2399      &  sccorname(:ilen(sccorname))
2400       write (iout,*) "Bond & inertia constant file    : ",
2401      &  bondname(:ilen(bondname))
2402       write (iout,*) "Bending-torsion parameter file  : ",
2403      &  thetname(:ilen(thetname))
2404       write (iout,*) "Rotamer parameter file          : ",
2405      &  rotname(:ilen(rotname))
2406       write (iout,*) "Thetpdb parameter file          : ",
2407      &  thetname_pdb(:ilen(thetname_pdb))
2408       write (iout,*) "Threading database              : ",
2409      &  patname(:ilen(patname))
2410       if (lentmp.ne.0)
2411      &write (iout,*)" DIRTMP                          : ",
2412      &  tmpdir(:lentmp)
2413       write (iout,'(80(1h-))')
2414       endif
2415       return
2416       end
2417 c----------------------------------------------------------------------------
2418       subroutine card_concat(card)
2419       implicit real*8 (a-h,o-z)
2420       include 'DIMENSIONS'
2421       include 'COMMON.IOUNITS'
2422       character*(*) card
2423       character*80 karta,ucase
2424       external ilen
2425       read (inp,'(a)') karta
2426       karta=ucase(karta)
2427       card=' '
2428       do while (karta(80:80).eq.'&')
2429         card=card(:ilen(card)+1)//karta(:79)
2430         read (inp,'(a)') karta
2431         karta=ucase(karta)
2432       enddo
2433       card=card(:ilen(card)+1)//karta
2434       return
2435       end
2436 c------------------------------------------------------------------------------
2437       subroutine readrst
2438       implicit none
2439       include 'DIMENSIONS'
2440       include 'COMMON.CHAIN'
2441       include 'COMMON.IOUNITS'
2442       include 'COMMON.CONTROL'
2443       include 'COMMON.MD'
2444       include 'COMMON.QRESTR'
2445 #ifdef FIVEDIAG
2446        include 'COMMON.LAGRANGE.5diag'
2447 #else
2448        include 'COMMON.LAGRANGE'
2449 #endif
2450       integer i,j
2451       open(irest2,file=rest2name,status='unknown')
2452       read(irest2,*) totT,EK,potE,totE,t_bath
2453       totTafm=totT
2454       do i=0,2*nres-1
2455          read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2456       enddo
2457       do i=0,2*nres-1
2458          read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2459       enddo
2460       if(usampl) then
2461              read (irest2,*) iset
2462       endif
2463       close(irest2)
2464       return
2465       end
2466 c------------------------------------------------------------------------------
2467       subroutine read_fragments
2468       implicit none
2469       include 'DIMENSIONS'
2470 #ifdef MPI
2471       include 'mpif.h'
2472 #endif
2473       include 'COMMON.SETUP'
2474       include 'COMMON.CHAIN'
2475       include 'COMMON.IOUNITS'
2476       include 'COMMON.MD'
2477       include 'COMMON.QRESTR'
2478       include 'COMMON.CONTROL'
2479       integer i
2480       read(inp,*) nset,nfrag,npair,nfrag_back
2481       loc_qlike=(nfrag_back.lt.0)
2482       nfrag_back=iabs(nfrag_back)
2483       if(me.eq.king.or..not.out1file)
2484      & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2485      &  " nfrag_back",nfrag_back," loc_qlike",loc_qlike
2486       do iset=1,nset
2487          read(inp,*) mset(iset)
2488        do i=1,nfrag
2489          read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2490      &     qinfrag(i,iset)
2491          if(me.eq.king.or..not.out1file)
2492      &    write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2493      &     ifrag(2,i,iset), qinfrag(i,iset)
2494        enddo
2495        do i=1,npair
2496         read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2497      &    qinpair(i,iset)
2498         if(me.eq.king.or..not.out1file)
2499      &   write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2500      &    ipair(2,i,iset), qinpair(i,iset)
2501        enddo
2502        if (loc_qlike) then
2503        do i=1,nfrag_back
2504         read(inp,*) wfrag_back(1,i,iset),qin_back(1,i,iset),
2505      &     wfrag_back(2,i,iset),qin_back(2,i,iset),
2506      &     wfrag_back(3,i,iset),qin_back(3,i,iset),
2507      &     ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2508         if(me.eq.king.or..not.out1file)
2509      &   write(iout,*) "A",i,wfrag_back(1,i,iset),qin_back(2,i,iset),
2510      &   wfrag_back(2,i,iset),qin_back(3,i,iset),
2511      &   wfrag_back(3,i,iset),qin_back(3,i,iset),
2512      &   ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2513        enddo
2514        else
2515        do i=1,nfrag_back
2516         read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2517      &     wfrag_back(3,i,iset),
2518      &     ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2519         if(me.eq.king.or..not.out1file)
2520      &   write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2521      &   wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2522        enddo
2523        endif
2524       enddo
2525       return
2526       end
2527 C---------------------------------------------------------------------------
2528       subroutine read_afminp
2529       implicit none
2530       include 'DIMENSIONS'
2531 #ifdef MPI
2532       include 'mpif.h'
2533 #endif
2534       include 'COMMON.SETUP'
2535       include 'COMMON.CONTROL'
2536       include 'COMMON.CHAIN'
2537       include 'COMMON.IOUNITS'
2538       include 'COMMON.SBRIDGE'
2539       character*320 afmcard
2540       integer i
2541 c      print *, "wchodze"
2542       call card_concat(afmcard)
2543       call readi(afmcard,"BEG",afmbeg,0)
2544       call readi(afmcard,"END",afmend,0)
2545       call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
2546       call reada(afmcard,"VEL",velAFMconst,0.0d0)
2547 c      print *,'FORCE=' ,forceAFMconst
2548 CCCC NOW PROPERTIES FOR AFM
2549        distafminit=0.0d0
2550        do i=1,3
2551         distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
2552        enddo
2553         distafminit=dsqrt(distafminit)
2554 c        print *,'initdist',distafminit
2555       return
2556       end
2557 c-------------------------------------------------------------------------------
2558       subroutine read_saxs_constr
2559       implicit none
2560       include 'DIMENSIONS'
2561 #ifdef MPI
2562       include 'mpif.h'
2563 #endif
2564       include 'COMMON.SETUP'
2565       include 'COMMON.CONTROL'
2566       include 'COMMON.SAXS'
2567       include 'COMMON.CHAIN'
2568       include 'COMMON.IOUNITS'
2569       include 'COMMON.SBRIDGE'
2570       double precision cm(3),cnorm
2571       integer i,j
2572 c      read(inp,*) nsaxs
2573       write (iout,*) "Calling read_saxs nsaxs",nsaxs
2574       call flush(iout)
2575       if (saxs_mode.eq.0) then
2576 c SAXS distance distribution
2577       do i=1,nsaxs
2578         read(inp,*) distsaxs(i),Psaxs(i)
2579       enddo
2580       Cnorm = 0.0d0
2581       do i=1,nsaxs
2582         Cnorm = Cnorm + Psaxs(i)
2583       enddo
2584       write (iout,*) "Cnorm",Cnorm
2585       do i=1,nsaxs
2586         Psaxs(i)=Psaxs(i)/Cnorm
2587       enddo
2588       write (iout,*) "Normalized distance distribution from SAXS"
2589       do i=1,nsaxs
2590         write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
2591       enddo
2592       Wsaxs0=0.0d0
2593       do i=1,nsaxs
2594         Wsaxs0=Wsaxs0-Psaxs(i)*dlog(Psaxs(i))
2595       enddo
2596       write (iout,*) "Wsaxs0",Wsaxs0
2597       else
2598 c SAXS "spheres".
2599       do i=1,nsaxs
2600         read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3)
2601       enddo
2602       do j=1,3
2603         cm(j)=0.0d0
2604       enddo
2605       do i=1,nsaxs
2606         do j=1,3
2607           cm(j)=cm(j)+Csaxs(j,i)
2608         enddo
2609       enddo
2610       do j=1,3
2611         cm(j)=cm(j)/nsaxs
2612       enddo
2613       do i=1,nsaxs
2614         do j=1,3
2615           Csaxs(j,i)=Csaxs(j,i)-cm(j)
2616         enddo
2617       enddo
2618       write (iout,*) "SAXS sphere coordinates"
2619       do i=1,nsaxs
2620         write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3)
2621       enddo
2622       endif
2623       return
2624       end
2625
2626 c-------------------------------------------------------------------------------
2627       subroutine read_dist_constr
2628       implicit none
2629       include 'DIMENSIONS'
2630 #ifdef MPI
2631       include 'mpif.h'
2632 #endif
2633       include 'COMMON.SETUP'
2634       include 'COMMON.CONTROL'
2635       include 'COMMON.CHAIN'
2636       include 'COMMON.IOUNITS'
2637       include 'COMMON.SBRIDGE'
2638       include 'COMMON.INTERACT'
2639       integer i,j,k,ii,jj,itemp,link_type,iiend,jjend,kk
2640       integer nfrag_,npair_,ndist_,ifrag_(2,100),ipair_(2,1000)
2641       double precision wfrag_(100),wpair_(1000)
2642       double precision ddjk,dist,dist_cut,fordepthmax
2643       character*5000 controlcard
2644       logical normalize,next
2645       integer restr_type
2646       double precision scal_bfac
2647       double precision xlink(4,0:4) /
2648 c           a          b       c     sigma
2649      &   0.0d0,0.0d0,0.0d0,0.0d0,                             ! default, no xlink potential
2650      &   0.00305218d0,9.46638d0,4.68901d0,4.74347d0,          ! ZL
2651      &   0.00214928d0,12.7517d0,0.00375009d0,6.13477d0,       ! ADH
2652      &   0.00184547d0,11.2678d0,0.00140292d0,7.00868d0,       ! PDH
2653      &   0.000161786d0,6.29273d0,4.40993d0,7.13956d0    /     ! DSS
2654 c      print *, "WCHODZE"
2655       write (iout,*) "Calling read_dist_constr"
2656 c      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2657 c      call flush(iout)
2658       restr_on_coord=.false.
2659       next=.true.
2660
2661       npeak=0
2662       ipeak=0
2663       nhpb_peak=0
2664
2665       DO WHILE (next)
2666
2667       call card_concat(controlcard)
2668       next = index(controlcard,"NEXT").gt.0
2669       call readi(controlcard,"RESTR_TYPE",restr_type,constr_dist)
2670       write (iout,*) "restr_type",restr_type
2671       call readi(controlcard,"NFRAG",nfrag_,0)
2672       call readi(controlcard,"NFRAG",nfrag_,0)
2673       call readi(controlcard,"NPAIR",npair_,0)
2674       call readi(controlcard,"NDIST",ndist_,0)
2675       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2676       call reada(controlcard,'SCAL_BFAC',scal_bfac,1.0d0)
2677       if (restr_type.eq.10)
2678      &  call reada(controlcard,'WBOLTZD',wboltzd,0.591d0)
2679       if (restr_type.eq.12)
2680      &  call reada(controlcard,'SCAL_PEAK',scal_peak,5.0d0)
2681       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2682       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2683       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2684       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2685       normalize = index(controlcard,"NORMALIZE").gt.0
2686       write (iout,*) "WBOLTZD",wboltzd
2687       write (iout,*) "SCAL_PEAK",scal_peak
2688       write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2689       write (iout,*) "IFRAG"
2690       do i=1,nfrag_
2691         write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2692       enddo
2693       write (iout,*) "IPAIR"
2694       do i=1,npair_
2695         write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2696       enddo
2697       if (nfrag_.gt.0 .or. restr_type.eq.4 .or. restr_type.eq.5)
2698      & write (iout,*)
2699      &   "Distance restraints as generated from reference structure"
2700       do i=1,nfrag_
2701         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2702         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2703      &    ifrag_(2,i)=nstart_sup+nsup-1
2704 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2705 c        call flush(iout)
2706         if (wfrag_(i).eq.0.0d0) cycle
2707         do j=ifrag_(1,i),ifrag_(2,i)-1
2708           do k=j+1,ifrag_(2,i)
2709 c            write (iout,*) "j",j," k",k
2710             ddjk=dist(j,k)
2711             if (restr_type.eq.1) then
2712               nhpb=nhpb+1
2713               irestr_type(nhpb)=1
2714               ihpb(nhpb)=j
2715               jhpb(nhpb)=k
2716               dhpb(nhpb)=ddjk
2717               forcon(nhpb)=wfrag_(i)
2718             else if (constr_dist.eq.2) then
2719               if (ddjk.le.dist_cut) then
2720                 nhpb=nhpb+1
2721                 irestr_type(nhpb)=1
2722                 ihpb(nhpb)=j
2723                 jhpb(nhpb)=k
2724                 dhpb(nhpb)=ddjk
2725                 forcon(nhpb)=wfrag_(i)
2726               endif
2727             else if (restr_type.eq.3) then
2728               nhpb=nhpb+1
2729               irestr_type(nhpb)=1
2730               ihpb(nhpb)=j
2731               jhpb(nhpb)=k
2732               dhpb(nhpb)=ddjk
2733               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2734             endif
2735 #ifdef MPI
2736             if (.not.out1file .or. me.eq.king)
2737      &      write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2738      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2739 #else
2740             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2741      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2742 #endif
2743           enddo
2744         enddo
2745       enddo
2746       do i=1,npair_
2747         if (wpair_(i).eq.0.0d0) cycle
2748         ii = ipair_(1,i)
2749         jj = ipair_(2,i)
2750         if (ii.gt.jj) then
2751           itemp=ii
2752           ii=jj
2753           jj=itemp
2754         endif
2755         do j=ifrag_(1,ii),ifrag_(2,ii)
2756           do k=ifrag_(1,jj),ifrag_(2,jj)
2757             ddjk=dist(j,k)
2758             if (restr_type.eq.1) then
2759               nhpb=nhpb+1
2760               irestr_type(nhpb)=1
2761               ihpb(nhpb)=j
2762               jhpb(nhpb)=k
2763               dhpb(nhpb)=ddjk
2764               forcon(nhpb)=wpair_(i)
2765             else if (constr_dist.eq.2) then
2766               if (ddjk.le.dist_cut) then
2767                 nhpb=nhpb+1
2768                 irestr_type(nhpb)=1
2769                 ihpb(nhpb)=j
2770                 jhpb(nhpb)=k
2771                 dhpb(nhpb)=ddjk
2772                 forcon(nhpb)=wpair_(i)
2773               endif
2774             else if (restr_type.eq.3) then
2775               nhpb=nhpb+1
2776               irestr_type(nhpb)=1
2777               ihpb(nhpb)=j
2778               jhpb(nhpb)=k
2779               dhpb(nhpb)=ddjk
2780               forcon(nhpb)=wpair_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2781             endif
2782 #ifdef MPI
2783             if (.not.out1file .or. me.eq.king)
2784      &      write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2785      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2786 #else
2787             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2788      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2789 #endif
2790           enddo
2791         enddo
2792       enddo
2793
2794 c      print *,ndist_
2795       write (iout,*) "Distance restraints as read from input"
2796       do i=1,ndist_
2797         if (restr_type.eq.12) then
2798           read (inp,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
2799      &    dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
2800      &    ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
2801      &    fordepth_peak(nhpb_peak+1),npeak
2802 c          write(iout,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
2803 c     &    dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
2804 c     &    ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
2805 c     &    fordepth_peak(nhpb_peak+1),npeak
2806           if (forcon_peak(nhpb_peak+1).le.0.0d0.or.
2807      &      fordepth_peak(nhpb_peak+1).le.0.0d0)cycle
2808           nhpb_peak=nhpb_peak+1
2809           irestr_type_peak(nhpb_peak)=12
2810           if (ipeak(1,npeak).eq.0) ipeak(1,npeak)=i
2811           ipeak(2,npeak)=i
2812 #ifdef MPI
2813           if (.not.out1file .or. me.eq.king)
2814      &    write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2815      &     nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
2816      &     ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
2817      &     dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
2818      &     fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
2819 #else
2820           write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2821      &     nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
2822      &     ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
2823      &     dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
2824      &     fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
2825 #endif
2826           if (ibecarb_peak(nhpb_peak).eq.3) then
2827             jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
2828           else if (ibecarb_peak(nhpb_peak).eq.2) then
2829             ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
2830           else if (ibecarb_peak(nhpb_peak).eq.1) then
2831             ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
2832             jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
2833           endif
2834         else if (restr_type.eq.11) then
2835           read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
2836      &     dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
2837 c        fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
2838           if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle
2839           nhpb=nhpb+1
2840           irestr_type(nhpb)=11
2841 #ifdef MPI
2842           if (.not.out1file .or. me.eq.king)
2843      &    write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2844      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2845      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
2846 #else
2847           write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2848      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2849      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
2850 #endif
2851 c          if (ibecarb(nhpb).gt.0) then
2852 c            ihpb(nhpb)=ihpb(nhpb)+nres
2853 c            jhpb(nhpb)=jhpb(nhpb)+nres
2854 c          endif
2855          if (ibecarb(nhpb).eq.3) then
2856             ihpb(nhpb)=ihpb(nhpb)+nres
2857           else if (ibecarb(nhpb).eq.2) then
2858             ihpb(nhpb)=ihpb(nhpb)+nres
2859           else if (ibecarb(nhpb).eq.1) then
2860             ihpb(nhpb)=ihpb(nhpb)+nres
2861             jhpb(nhpb)=jhpb(nhpb)+nres
2862           endif
2863         else if (restr_type.eq.10) then
2864 c Cross-lonk Markov-like potential
2865           call card_concat(controlcard)
2866           call readi(controlcard,"ILINK",ihpb(nhpb+1),0)
2867           call readi(controlcard,"JLINK",jhpb(nhpb+1),0)
2868           ibecarb(nhpb+1)=0
2869           if (index(controlcard,"BETA").gt.0) ibecarb(nhpb+1)=1
2870           if (ihpb(nhpb+1).eq.0 .or. jhpb(nhpb+1).eq.0) cycle
2871           if (index(controlcard,"ZL").gt.0) then
2872             link_type=1
2873           else if (index(controlcard,"ADH").gt.0) then
2874             link_type=2
2875           else if (index(controlcard,"PDH").gt.0) then
2876             link_type=3
2877           else if (index(controlcard,"DSS").gt.0) then
2878             link_type=4
2879           else
2880             link_type=0
2881           endif
2882           call reada(controlcard,"AXLINK",dhpb(nhpb+1),
2883      &       xlink(1,link_type))
2884           call reada(controlcard,"BXLINK",dhpb1(nhpb+1),
2885      &       xlink(2,link_type))
2886           call reada(controlcard,"CXLINK",fordepth(nhpb+1),
2887      &       xlink(3,link_type))
2888           call reada(controlcard,"SIGMA",forcon(nhpb+1),
2889      &       xlink(4,link_type))
2890           call reada(controlcard,"SCORE",xlscore(nhpb+1),1.0d0)
2891 c          read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),ibecarb(nhpb+1),
2892 c     &      dhpb(nhpb+1),dhpb1(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
2893           if (forcon(nhpb+1).le.0.0d0 .or.
2894      &       (dhpb(nhpb+1).eq.0 .and. dhpb1(nhpb+1).eq.0)) cycle
2895           nhpb=nhpb+1
2896           irestr_type(nhpb)=10
2897           if (ibecarb(nhpb).eq.3) then
2898             jhpb(nhpb)=jhpb(nhpb)+nres
2899           else if (ibecarb(nhpb).eq.2) then
2900             ihpb(nhpb)=ihpb(nhpb)+nres
2901           else if (ibecarb(nhpb).eq.1) then
2902             ihpb(nhpb)=ihpb(nhpb)+nres
2903             jhpb(nhpb)=jhpb(nhpb)+nres
2904           endif
2905 #ifdef MPI
2906           if (.not.out1file .or. me.eq.king)
2907      &    write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2908      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2909      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2910      &     irestr_type(nhpb)
2911 #else
2912           write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2913      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2914      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2915      &     irestr_type(nhpb)
2916 #endif
2917         else
2918 C        print *,"in else"
2919           read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
2920      &     dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1)
2921           if (forcon(nhpb+1).gt.0.0d0) then
2922           nhpb=nhpb+1
2923           if (dhpb1(nhpb).eq.0.0d0) then
2924             irestr_type(nhpb)=1
2925           else
2926             irestr_type(nhpb)=2
2927           endif
2928           if (ibecarb(nhpb).eq.3) then
2929             jhpb(nhpb)=jhpb(nhpb)+nres
2930           else if (ibecarb(nhpb).eq.2) then
2931             ihpb(nhpb)=ihpb(nhpb)+nres
2932           else if (ibecarb(nhpb).eq.1) then
2933             ihpb(nhpb)=ihpb(nhpb)+nres
2934             jhpb(nhpb)=jhpb(nhpb)+nres
2935           endif
2936           if (dhpb(nhpb).eq.0.0d0)
2937      &       dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2938         endif
2939 #ifdef MPI
2940           if (.not.out1file .or. me.eq.king)
2941      &    write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
2942      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
2943 #else
2944           write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
2945      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
2946 #endif
2947         endif
2948 C        read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2949 C        if (forcon(nhpb+1).gt.0.0d0) then
2950 C          nhpb=nhpb+1
2951 C          dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2952       enddo
2953
2954       if (restr_type.eq.4) then
2955         write (iout,*) "The BFAC array"
2956         do i=nnt,nct
2957           write (iout,'(i5,f10.5)') i,bfac(i)
2958         enddo
2959         do i=nnt,nct
2960           if (itype(i).eq.ntyp1) cycle
2961           do j=nnt,i-1
2962             if (itype(j).eq.ntyp1) cycle
2963             if (itype(i).eq.10) then
2964               iiend=0
2965             else
2966               iiend=1
2967             endif
2968             if (itype(j).eq.10) then
2969               jjend=0
2970             else
2971               jjend=1
2972             endif
2973             kk=0
2974             do ii=0,iiend
2975             do jj=0,jjend
2976             nhpb=nhpb+1
2977             irestr_type(nhpb)=1
2978             forcon(nhpb)=scal_bfac**2/(bfac(i)**2+bfac(j)**2)
2979             irestr_type(nhpb)=1
2980             ibecarb(nhpb)=kk
2981             if (ibecarb(nhpb).gt.0) ibecarb(nhpb)=4-ibecarb(nhpb)
2982             ihpb(nhpb)=i+nres*ii
2983             jhpb(nhpb)=j+nres*jj
2984             dhpb(nhpb)=dist(i+nres*ii,j+nres*jj)
2985 #ifdef MPI
2986             if (.not.out1file .or. me.eq.king) then
2987             write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2988      &       nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2989      &       dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2990      &       irestr_type(nhpb)
2991             endif
2992 #else
2993             write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2994      &       nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2995      &       dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2996      &       irestr_type(nhpb)
2997 #endif
2998             kk=kk+1
2999           enddo
3000           enddo
3001           enddo
3002         enddo
3003       endif
3004
3005       if (restr_type.eq.5) then
3006         restr_on_coord=.true.
3007         do i=nnt,nct
3008           if (itype(i).eq.ntyp1) cycle
3009           bfac(i)=(scal_bfac/bfac(i))**2
3010         enddo
3011       endif
3012
3013       ENDDO ! next
3014
3015       fordepthmax=0.0d0
3016       if (normalize) then
3017         do i=nss+1,nhpb
3018           if (irestr_type(i).eq.11.and.fordepth(i).gt.fordepthmax)
3019      &      fordepthmax=fordepth(i)
3020         enddo
3021         do i=nss+1,nhpb
3022           if (irestr_type(i).eq.11) fordepth(i)=fordepth(i)/fordepthmax
3023         enddo
3024       endif
3025       return
3026       end
3027 c-------------------------------------------------------------------------------
3028       subroutine read_constr_homology
3029       implicit none
3030       include 'DIMENSIONS'
3031 #ifdef MPI
3032       include 'mpif.h'
3033 #endif
3034       include 'COMMON.SETUP'
3035       include 'COMMON.CONTROL'
3036       include 'COMMON.HOMOLOGY'
3037       include 'COMMON.CHAIN'
3038       include 'COMMON.IOUNITS'
3039       include 'COMMON.MD'
3040       include 'COMMON.QRESTR'
3041       include 'COMMON.GEO'
3042       include 'COMMON.INTERACT'
3043       include 'COMMON.NAMES'
3044 c
3045 c For new homol impl
3046 c
3047       include 'COMMON.VAR'
3048 c
3049
3050 c     double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
3051 c    &                 dist_cut
3052 c     common /przechowalnia/ odl_temp(maxres,maxres,max_template),
3053 c    &    sigma_odl_temp(maxres,maxres,max_template)
3054       character*2 kic2
3055       character*24 model_ki_dist, model_ki_angle
3056       character*500 controlcard
3057       integer ki,i,ii,j,k,l,ii_in_use(maxdim),i_tmp,idomain_tmp,irec,
3058      & ik,iistart,nres_temp
3059       integer ilen
3060       external ilen
3061       logical liiflag,lfirst
3062       integer i01,i10
3063 c
3064 c     FP - Nov. 2014 Temporary specifications for new vars
3065 c
3066       double precision rescore_tmp,x12,y12,z12,rescore2_tmp,
3067      &                 rescore3_tmp
3068       double precision, dimension (max_template,maxres) :: rescore
3069       double precision, dimension (max_template,maxres) :: rescore2
3070       double precision, dimension (max_template,maxres) :: rescore3
3071       double precision distal
3072       character*24 tpl_k_rescore
3073       character*256 pdbfile
3074 c -----------------------------------------------------------------
3075 c Reading multiple PDB ref structures and calculation of retraints
3076 c not using pre-computed ones stored in files model_ki_{dist,angle}
3077 c FP (Nov., 2014)
3078 c -----------------------------------------------------------------
3079 c
3080 c
3081 c Alternative: reading from input
3082       call card_concat(controlcard)
3083       call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
3084       call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
3085       call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
3086       call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
3087       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
3088       call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
3089       call readi(controlcard,"HOMOL_NSET",homol_nset,1)
3090       read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
3091       start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
3092       if(.not.read2sigma.and.start_from_model) then
3093           if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
3094      &      write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA'
3095           start_from_model=.false.
3096           iranconf=(indpdb.le.0)
3097       endif
3098       if(start_from_model .and. (me.eq.king .or. .not. out1file))
3099      &    write(iout,*) 'START_FROM_MODELS is ON'
3100 c      if(start_from_model .and. rest) then
3101 c        if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3102 c         write(iout,*) 'START_FROM_MODELS is OFF'
3103 c         write(iout,*) 'remove restart keyword from input'
3104 c        endif
3105 c      endif
3106       if (start_from_model) nmodel_start=constr_homology
3107       if (homol_nset.gt.1)then
3108          call card_concat(controlcard)
3109          read(controlcard,*) (waga_homology(i),i=1,homol_nset)
3110          if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3111           write(iout,*) "iset homology_weight "
3112           do i=1,homol_nset
3113            write(iout,*) i,waga_homology(i)
3114           enddo
3115          endif
3116          iset=mod(kolor,homol_nset)+1
3117       else
3118        iset=1
3119        waga_homology(1)=1.0
3120       endif
3121
3122 cd      write (iout,*) "nnt",nnt," nct",nct
3123 cd      call flush(iout)
3124
3125
3126       lim_odl=0
3127       lim_dih=0
3128 c
3129 c      write(iout,*) 'nnt=',nnt,'nct=',nct
3130 c
3131       do i = nnt,nct
3132         do k=1,constr_homology
3133          idomain(k,i)=0
3134         enddo
3135       enddo
3136
3137       ii=0
3138       do i = nnt,nct-2
3139         do j=i+2,nct
3140         ii=ii+1
3141         ii_in_use(ii)=0
3142         enddo
3143       enddo
3144
3145       if (read_homol_frag) then
3146         call read_klapaucjusz
3147       else
3148
3149       do k=1,constr_homology
3150
3151         read(inp,'(a)') pdbfile
3152         pdbfiles_chomo(k)=pdbfile
3153         if(me.eq.king .or. .not. out1file)
3154      &    write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
3155      &  pdbfile(:ilen(pdbfile))
3156         open(ipdbin,file=pdbfile,status='old',err=33)
3157         goto 34
3158   33    write (iout,'(a,5x,a)') 'Error opening PDB file',
3159      &  pdbfile(:ilen(pdbfile))
3160         stop
3161   34    continue
3162 c        print *,'Begin reading pdb data'
3163 c
3164 c Files containing res sim or local scores (former containing sigmas)
3165 c
3166
3167         write(kic2,'(bz,i2.2)') k
3168
3169         tpl_k_rescore="template"//kic2//".sco"
3170
3171         unres_pdb=.false.
3172         nres_temp=nres
3173         if (read2sigma) then
3174           call readpdb_template(k)
3175         else
3176           call readpdb
3177         endif
3178         nres_chomo(k)=nres
3179         nres=nres_temp
3180 c
3181 c     Distance restraints
3182 c
3183 c          ... --> odl(k,ii)
3184 C Copy the coordinates from reference coordinates (?)
3185         do i=1,2*nres_chomo(k)
3186           do j=1,3
3187             c(j,i)=cref(j,i)
3188 c           write (iout,*) "c(",j,i,") =",c(j,i)
3189           enddo
3190         enddo
3191 c
3192 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
3193 c
3194 c         write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
3195           open (ientin,file=tpl_k_rescore,status='old')
3196           if (nnt.gt.1) rescore(k,1)=0.0d0
3197           do irec=nnt,nct ! loop for reading res sim
3198             if (read2sigma) then
3199              read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
3200      &                                rescore3_tmp,idomain_tmp
3201              i_tmp=i_tmp+nnt-1
3202              idomain(k,i_tmp)=idomain_tmp
3203              rescore(k,i_tmp)=rescore_tmp
3204              rescore2(k,i_tmp)=rescore2_tmp
3205              rescore3(k,i_tmp)=rescore3_tmp
3206              if (.not. out1file .or. me.eq.king)
3207      &        write(iout,'(a7,i5,3f10.5,i5)') "rescore",
3208      &                      i_tmp,rescore2_tmp,rescore_tmp,
3209      &                                rescore3_tmp,idomain_tmp
3210             else
3211              idomain(k,irec)=1
3212              read (ientin,*,end=1401) rescore_tmp
3213
3214 c           rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
3215              rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
3216 c           write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
3217             endif
3218           enddo
3219  1401   continue
3220         close (ientin)
3221         if (waga_dist.ne.0.0d0) then
3222           ii=0
3223           do i = nnt,nct-2
3224             do j=i+2,nct
3225
3226               x12=c(1,i)-c(1,j)
3227               y12=c(2,i)-c(2,j)
3228               z12=c(3,i)-c(3,j)
3229               distal=dsqrt(x12*x12+y12*y12+z12*z12)
3230 c              write (iout,*) k,i,j,distal,dist2_cut
3231
3232             if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3233      &            .and. distal.le.dist2_cut ) then
3234
3235               ii=ii+1
3236               ii_in_use(ii)=1
3237               l_homo(k,ii)=.true.
3238
3239 c             write (iout,*) "k",k
3240 c             write (iout,*) "i",i," j",j," constr_homology",
3241 c    &                       constr_homology
3242               ires_homo(ii)=i
3243               jres_homo(ii)=j
3244               odl(k,ii)=distal
3245               if (read2sigma) then
3246                 sigma_odl(k,ii)=0
3247                 do ik=i,j
3248                  sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
3249                 enddo
3250                 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
3251                 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
3252      &        sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3253               else
3254                 if (odl(k,ii).le.dist_cut) then
3255                  sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
3256                 else
3257 #ifdef OLDSIGMA
3258                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3259      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
3260 #else
3261                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3262      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3263 #endif
3264                 endif
3265               endif
3266               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
3267             else
3268               ii=ii+1
3269               l_homo(k,ii)=.false.
3270             endif
3271             enddo
3272           enddo
3273         lim_odl=ii
3274         endif
3275 c        write (iout,*) "Distance restraints set"
3276 c        call flush(iout)
3277 c
3278 c     Theta, dihedral and SC retraints
3279 c
3280         if (waga_angle.gt.0.0d0) then
3281 c         open (ientin,file=tpl_k_sigma_dih,status='old')
3282 c         do irec=1,maxres-3 ! loop for reading sigma_dih
3283 c            read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
3284 c            if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
3285 c            sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
3286 c    &                            sigma_dih(k,i+nnt-1)
3287 c         enddo
3288 c1402   continue
3289 c         close (ientin)
3290           do i = nnt+3,nct
3291             if (idomain(k,i).eq.0) then
3292                sigma_dih(k,i)=0.0
3293                cycle
3294             endif
3295             dih(k,i)=phiref(i) ! right?
3296 c           read (ientin,*) sigma_dih(k,i) ! original variant
3297 c             write (iout,*) "dih(",k,i,") =",dih(k,i)
3298 c             write(iout,*) "rescore(",k,i,") =",rescore(k,i),
3299 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
3300 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2),
3301 c    &                      "rescore(",k,i-3,") =",rescore(k,i-3)
3302
3303             sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
3304      &                     rescore(k,i-2)+rescore(k,i-3))/4.0
3305 c            if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
3306 c           write (iout,*) "Raw sigmas for dihedral angle restraints"
3307 c           write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
3308 c           sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
3309 c                          rescore(k,i-2)*rescore(k,i-3)  !  right expression ?
3310 c   Instead of res sim other local measure of b/b str reliability possible
3311             if (sigma_dih(k,i).ne.0)
3312      &       sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
3313 c           sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
3314           enddo
3315           lim_dih=nct-nnt-2
3316         endif
3317 c        write (iout,*) "Dihedral angle restraints set"
3318 c        call flush(iout)
3319
3320         if (waga_theta.gt.0.0d0) then
3321 c         open (ientin,file=tpl_k_sigma_theta,status='old')
3322 c         do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
3323 c            read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
3324 c            sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
3325 c    &                              sigma_theta(k,i+nnt-1)
3326 c         enddo
3327 c1403   continue
3328 c         close (ientin)
3329
3330           do i = nnt+2,nct ! right? without parallel.
3331 c         do i = i=1,nres ! alternative for bounds acc to readpdb?
3332 c         do i=ithet_start,ithet_end ! with FG parallel.
3333              if (idomain(k,i).eq.0) then
3334               sigma_theta(k,i)=0.0
3335               cycle
3336              endif
3337              thetatpl(k,i)=thetaref(i)
3338 c            write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
3339 c            write(iout,*)  "rescore(",k,i,") =",rescore(k,i),
3340 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
3341 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2)
3342 c            read (ientin,*) sigma_theta(k,i) ! 1st variant
3343              sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
3344      &                        rescore(k,i-2))/3.0
3345 c             if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
3346              if (sigma_theta(k,i).ne.0)
3347      &       sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
3348
3349 c            sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
3350 c                             rescore(k,i-2) !  right expression ?
3351 c            sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
3352           enddo
3353         endif
3354 c        write (iout,*) "Angle restraints set"
3355 c        call flush(iout)
3356
3357         if (waga_d.gt.0.0d0) then
3358 c       open (ientin,file=tpl_k_sigma_d,status='old')
3359 c         do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
3360 c            read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
3361 c            sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
3362 c    &                          sigma_d(k,i+nnt-1)
3363 c         enddo
3364 c1404   continue
3365
3366           do i = nnt,nct ! right? without parallel.
3367 c         do i=2,nres-1 ! alternative for bounds acc to readpdb?
3368 c         do i=loc_start,loc_end ! with FG parallel.
3369                if (itype(i).eq.10) cycle
3370                if (idomain(k,i).eq.0 ) then
3371                   sigma_d(k,i)=0.0
3372                   cycle
3373                endif
3374                xxtpl(k,i)=xxref(i)
3375                yytpl(k,i)=yyref(i)
3376                zztpl(k,i)=zzref(i)
3377 c              write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
3378 c              write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
3379 c              write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
3380 c              write(iout,*)  "rescore(",k,i,") =",rescore(k,i)
3381                sigma_d(k,i)=rescore3(k,i) !  right expression ?
3382                if (sigma_d(k,i).ne.0)
3383      &          sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
3384
3385 c              sigma_d(k,i)=hmscore(k)*rescore(k,i) !  right expression ?
3386 c              sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
3387 c              read (ientin,*) sigma_d(k,i) ! 1st variant
3388           enddo
3389         endif
3390       enddo
3391 c      write (iout,*) "SC restraints set"
3392 c      call flush(iout)
3393 c
3394 c remove distance restraints not used in any model from the list
3395 c shift data in all arrays
3396 c
3397 c      write (iout,*) "waga_dist",waga_dist," nnt",nnt," nct",nct
3398       if (waga_dist.ne.0.0d0) then
3399         ii=0
3400         liiflag=.true.
3401         lfirst=.true.
3402         do i=nnt,nct-2
3403          do j=i+2,nct
3404           ii=ii+1
3405 c          if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3406 c     &            .and. distal.le.dist2_cut ) then
3407 c          write (iout,*) "i",i," j",j," ii",ii
3408 c          call flush(iout)
3409           if (ii_in_use(ii).eq.0.and.liiflag.or.
3410      &     ii_in_use(ii).eq.1.and.liiflag.and.ii.eq.lim_odl) then
3411             liiflag=.false.
3412             i10=ii
3413             if (lfirst) then
3414               lfirst=.false.
3415               iistart=ii
3416             else
3417               if(i10.eq.lim_odl) i10=i10+1
3418               do ki=0,i10-i01-1
3419                ires_homo(iistart+ki)=ires_homo(ki+i01)
3420                jres_homo(iistart+ki)=jres_homo(ki+i01)
3421                ii_in_use(iistart+ki)=ii_in_use(ki+i01)
3422                do k=1,constr_homology
3423                 odl(k,iistart+ki)=odl(k,ki+i01)
3424                 sigma_odl(k,iistart+ki)=sigma_odl(k,ki+i01)
3425                 l_homo(k,iistart+ki)=l_homo(k,ki+i01)
3426                enddo
3427               enddo
3428               iistart=iistart+i10-i01
3429             endif
3430           endif
3431           if (ii_in_use(ii).ne.0.and..not.liiflag) then
3432              i01=ii
3433              liiflag=.true.
3434           endif
3435          enddo
3436         enddo
3437         lim_odl=iistart-1
3438       endif
3439 c      write (iout,*) "Removing distances completed"
3440 c      call flush(iout)
3441       endif ! .not. klapaucjusz
3442
3443       if (constr_homology.gt.0) call homology_partition
3444 c      write (iout,*) "After homology_partition"
3445 c      call flush(iout)
3446       if (constr_homology.gt.0) call init_int_table
3447 c      write (iout,*) "After init_int_table"
3448 c      call flush(iout)
3449 c      write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
3450 c      write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
3451 c
3452 c Print restraints
3453 c
3454       if (.not.out_template_restr) return
3455 cd      write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
3456       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3457        write (iout,*) "Distance restraints from templates"
3458        do ii=1,lim_odl
3459        write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
3460      &  ii,ires_homo(ii),jres_homo(ii),
3461      &  (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
3462      &  ki=1,constr_homology)
3463        enddo
3464        write (iout,*) "Dihedral angle restraints from templates"
3465        do i=nnt+3,nct
3466         write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3467      &      (rad2deg*dih(ki,i),
3468      &      rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
3469        enddo
3470        write (iout,*) "Virtual-bond angle restraints from templates"
3471        do i=nnt+2,nct
3472         write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3473      &      (rad2deg*thetatpl(ki,i),
3474      &      rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
3475        enddo
3476        write (iout,*) "SC restraints from templates"
3477        do i=nnt,nct
3478         write(iout,'(i5,100(4f8.2,4x))') i,
3479      &  (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
3480      &   1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
3481        enddo
3482       endif
3483 c -----------------------------------------------------------------
3484       return
3485       end
3486 c----------------------------------------------------------------------
3487 #ifdef WINIFL
3488       subroutine flush(iu)
3489       return
3490       end
3491 #endif
3492 #ifdef AIX
3493       subroutine flush(iu)
3494       call flush_(iu)
3495       return
3496       end
3497 #endif
3498 c------------------------------------------------------------------------------
3499       subroutine copy_to_tmp(source)
3500       implicit none
3501       include "DIMENSIONS"
3502       include "COMMON.IOUNITS"
3503       character*(*) source
3504       character* 256 tmpfile
3505       integer ilen
3506       external ilen
3507       logical ex
3508       tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
3509       inquire(file=tmpfile,exist=ex)
3510       if (ex) then
3511         write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
3512      &   " to temporary directory..."
3513         write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
3514         call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
3515       endif
3516       return
3517       end
3518 c------------------------------------------------------------------------------
3519       subroutine move_from_tmp(source)
3520       implicit none
3521       include "DIMENSIONS"
3522       include "COMMON.IOUNITS"
3523       character*(*) source
3524       integer ilen
3525       external ilen
3526       write (*,*) "Moving ",source(:ilen(source)),
3527      & " from temporary directory to working directory"
3528       write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
3529       call system("/bin/mv "//source(:ilen(source))//" "//curdir)
3530       return
3531       end
3532 c------------------------------------------------------------------------------
3533       subroutine random_init(seed)
3534 C
3535 C Initialize random number generator
3536 C
3537       implicit none
3538       include 'DIMENSIONS'
3539 #ifdef MPI
3540       include 'mpif.h'
3541       logical OKRandom, prng_restart
3542       real*8  r1
3543       integer iseed_array(4)
3544       integer error_msg,ierr
3545 #endif
3546       include 'COMMON.IOUNITS'
3547       include 'COMMON.TIME1'
3548       include 'COMMON.THREAD'
3549       include 'COMMON.SBRIDGE'
3550       include 'COMMON.CONTROL'
3551       include 'COMMON.MCM'
3552       include 'COMMON.MAP'
3553       include 'COMMON.HEADER'
3554       include 'COMMON.CSA'
3555       include 'COMMON.CHAIN'
3556       include 'COMMON.MUCA'
3557       include 'COMMON.MD'
3558       include 'COMMON.FFIELD'
3559       include 'COMMON.SETUP'
3560       integer i,iseed
3561       double precision seed,ran_number
3562       iseed=-dint(dabs(seed))
3563       if (iseed.eq.0) then
3564         write (iout,'(/80(1h*)/20x,a/80(1h*))')
3565      &    'Random seed undefined. The program will stop.'
3566         write (*,'(/80(1h*)/20x,a/80(1h*))')
3567      &    'Random seed undefined. The program will stop.'
3568 #ifdef MPI
3569         call mpi_finalize(mpi_comm_world,ierr)
3570 #endif
3571         stop 'Bad random seed.'
3572       endif
3573 #ifdef MPI
3574       if (fg_rank.eq.0) then
3575       seed=seed*(me+1)+1
3576 #ifdef AMD64
3577       if(me.eq.king)
3578      &  write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
3579       OKRandom = prng_restart(me,iseed)
3580 #else
3581       do i=1,4
3582        tmp=65536.0d0**(4-i)
3583        iseed_array(i) = dint(seed/tmp)
3584        seed=seed-iseed_array(i)*tmp
3585       enddo
3586       if(me.eq.king)
3587      & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
3588      &                 (iseed_array(i),i=1,4)
3589       write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
3590      &                 (iseed_array(i),i=1,4)
3591       OKRandom = prng_restart(me,iseed_array)
3592 #endif
3593       if (OKRandom) then
3594 c        r1 = prng_next(me)
3595         r1=ran_number(0.0D0,1.0D0)
3596         if(me.eq.king)
3597      &   write (iout,*) 'ran_num',r1
3598         if (r1.lt.0.0d0) OKRandom=.false.
3599       endif
3600       if (.not.OKRandom) then
3601         write (iout,*) 'PRNG IS NOT WORKING!!!'
3602         print *,'PRNG IS NOT WORKING!!!'
3603         if (me.eq.0) then
3604          call flush(iout)
3605          call mpi_abort(mpi_comm_world,error_msg,ierr)
3606          stop
3607         else
3608          write (iout,*) 'too many processors for parallel prng'
3609          write (*,*) 'too many processors for parallel prng'
3610          call flush(iout)
3611          stop
3612         endif
3613       endif
3614       endif
3615 #else
3616       call vrndst(iseed)
3617       write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
3618 #endif
3619       return
3620       end
3621 c----------------------------------------------------------------------
3622       subroutine read_klapaucjusz
3623       implicit none
3624       include 'DIMENSIONS'
3625 #ifdef MPI
3626       include 'mpif.h'
3627 #endif
3628       include 'COMMON.SETUP'
3629       include 'COMMON.CONTROL'
3630       include 'COMMON.HOMOLOGY'
3631       include 'COMMON.CHAIN'
3632       include 'COMMON.IOUNITS'
3633       include 'COMMON.MD'
3634       include 'COMMON.GEO'
3635       include 'COMMON.INTERACT'
3636       include 'COMMON.NAMES'
3637       character*256 fragfile
3638       integer ninclust(maxclust),inclust(max_template,maxclust),
3639      &  nresclust(maxclust),iresclust(maxres,maxclust),nclust
3640
3641       character*2 kic2
3642       character*24 model_ki_dist, model_ki_angle
3643       character*500 controlcard
3644       integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp,
3645      & ik,ll,ii,kk,iistart,iishift,lim_xx
3646       double precision distal
3647       logical lprn /.true./
3648       integer nres_temp
3649       integer ilen
3650       external ilen
3651       logical liiflag
3652 c
3653 c
3654       double precision rescore_tmp,x12,y12,z12,rescore2_tmp
3655       double precision, dimension (max_template,maxres) :: rescore
3656       double precision, dimension (max_template,maxres) :: rescore2
3657       character*24 tpl_k_rescore
3658       character*256 pdbfile
3659
3660 c
3661 c For new homol impl
3662 c
3663       include 'COMMON.VAR'
3664 c
3665       call getenv("FRAGFILE",fragfile)
3666       open(ientin,file=fragfile,status="old",err=10)
3667       read(ientin,*) constr_homology,nclust
3668       nmodel_start=constr_homology
3669       l_homo = .false.
3670       sigma_theta=0.0
3671       sigma_d=0.0
3672       sigma_dih=0.0
3673 c Read pdb files
3674       do k=1,constr_homology
3675         read(ientin,'(a)') pdbfile
3676         write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file',
3677      &  pdbfile(:ilen(pdbfile))
3678         pdbfiles_chomo(k)=pdbfile
3679         open(ipdbin,file=pdbfile,status='old',err=33)
3680         goto 34
3681   33    write (iout,'(a,5x,a)') 'Error opening PDB file',
3682      &  pdbfile(:ilen(pdbfile))
3683         stop
3684   34    continue
3685         unres_pdb=.false.
3686         nres_temp=nres
3687         call readpdb_template(k)
3688         nres_chomo(k)=nres
3689         nres=nres_temp
3690         do i=1,nres
3691           rescore(k,i)=0.2d0
3692           rescore2(k,i)=1.0d0
3693         enddo
3694       enddo
3695 c Read clusters
3696       do i=1,nclust
3697         read(ientin,*) ninclust(i),nresclust(i)
3698         read(ientin,*) (inclust(k,i),k=1,ninclust(i))
3699         read(ientin,*) (iresclust(k,i),k=1,nresclust(i))
3700       enddo
3701 c
3702 c Loop over clusters
3703 c
3704       do l=1,nclust
3705         do ll = 1,ninclust(l)
3706
3707         k = inclust(ll,l)
3708         do i=1,nres
3709           idomain(k,i)=0
3710         enddo
3711         do i=1,nresclust(l)
3712           if (nnt.gt.1)  then
3713             idomain(k,iresclust(i,l)+1) = 1
3714           else
3715             idomain(k,iresclust(i,l)) = 1
3716           endif
3717         enddo
3718 c
3719 c     Distance restraints
3720 c
3721 c          ... --> odl(k,ii)
3722 C Copy the coordinates from reference coordinates (?)
3723         nres_temp=nres
3724         nres=nres_chomo(k)
3725         do i=1,2*nres
3726           do j=1,3
3727             c(j,i)=chomo(j,i,k)
3728 c           write (iout,*) "c(",j,i,") =",c(j,i)
3729           enddo
3730         enddo
3731         call int_from_cart(.true.,.false.)
3732         call sc_loc_geom(.false.)
3733         do i=1,nres
3734           thetaref(i)=theta(i)
3735           phiref(i)=phi(i)
3736         enddo
3737         nres=nres_temp
3738         if (waga_dist.ne.0.0d0) then
3739           ii=0
3740           do i = nnt,nct-2
3741             do j=i+2,nct
3742
3743               x12=c(1,i)-c(1,j)
3744               y12=c(2,i)-c(2,j)
3745               z12=c(3,i)-c(3,j)
3746               distal=dsqrt(x12*x12+y12*y12+z12*z12)
3747 c              write (iout,*) k,i,j,distal,dist2_cut
3748
3749             if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3750      &            .and. distal.le.dist2_cut ) then
3751
3752               ii=ii+1
3753               ii_in_use(ii)=1
3754               l_homo(k,ii)=.true.
3755
3756 c             write (iout,*) "k",k
3757 c             write (iout,*) "i",i," j",j," constr_homology",
3758 c    &                       constr_homology
3759               ires_homo(ii)=i
3760               jres_homo(ii)=j
3761               odl(k,ii)=distal
3762               if (read2sigma) then
3763                 sigma_odl(k,ii)=0
3764                 do ik=i,j
3765                  sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
3766                 enddo
3767                 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
3768                 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
3769      &        sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3770               else
3771                 if (odl(k,ii).le.dist_cut) then
3772                  sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
3773                 else
3774 #ifdef OLDSIGMA
3775                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3776      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
3777 #else
3778                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3779      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3780 #endif
3781                 endif
3782               endif
3783               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
3784             else
3785               ii=ii+1
3786 c              l_homo(k,ii)=.false.
3787             endif
3788             enddo
3789           enddo
3790         lim_odl=ii
3791         endif
3792 c
3793 c     Theta, dihedral and SC retraints
3794 c
3795         if (waga_angle.gt.0.0d0) then
3796           do i = nnt+3,nct
3797             if (idomain(k,i).eq.0) then
3798 c               sigma_dih(k,i)=0.0
3799                cycle
3800             endif
3801             dih(k,i)=phiref(i)
3802             sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
3803      &                     rescore(k,i-2)+rescore(k,i-3))/4.0
3804 c            write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
3805 c     &       " sigma_dihed",sigma_dih(k,i)
3806             if (sigma_dih(k,i).ne.0)
3807      &       sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
3808           enddo
3809           lim_dih=nct-nnt-2
3810         endif
3811
3812         if (waga_theta.gt.0.0d0) then
3813           do i = nnt+2,nct
3814              if (idomain(k,i).eq.0) then
3815 c              sigma_theta(k,i)=0.0
3816               cycle
3817              endif
3818              thetatpl(k,i)=thetaref(i)
3819              sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
3820      &                        rescore(k,i-2))/3.0
3821              if (sigma_theta(k,i).ne.0)
3822      &       sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
3823           enddo
3824         endif
3825
3826         if (waga_d.gt.0.0d0) then
3827           do i = nnt,nct
3828                if (itype(i).eq.10) cycle
3829                if (idomain(k,i).eq.0 ) then
3830 c                  sigma_d(k,i)=0.0
3831                   cycle
3832                endif
3833                xxtpl(k,i)=xxref(i)
3834                yytpl(k,i)=yyref(i)
3835                zztpl(k,i)=zzref(i)
3836                sigma_d(k,i)=rescore(k,i)
3837                if (sigma_d(k,i).ne.0)
3838      &          sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
3839                if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
3840           enddo
3841         endif
3842       enddo ! l
3843       enddo ! ll
3844 c
3845 c remove distance restraints not used in any model from the list
3846 c shift data in all arrays
3847 c
3848       if (waga_dist.ne.0.0d0) then
3849         ii=0
3850         liiflag=.true.
3851         do i=nnt,nct-2
3852          do j=i+2,nct
3853           ii=ii+1
3854           if (ii_in_use(ii).eq.0.and.liiflag) then
3855             liiflag=.false.
3856             iistart=ii
3857           endif
3858           if (ii_in_use(ii).ne.0.and..not.liiflag.or.
3859      &                   .not.liiflag.and.ii.eq.lim_odl) then
3860              if (ii.eq.lim_odl) then
3861               iishift=ii-iistart+1
3862              else
3863               iishift=ii-iistart
3864              endif
3865              liiflag=.true.
3866              do ki=iistart,lim_odl-iishift
3867               ires_homo(ki)=ires_homo(ki+iishift)
3868               jres_homo(ki)=jres_homo(ki+iishift)
3869               ii_in_use(ki)=ii_in_use(ki+iishift)
3870               do k=1,constr_homology
3871                odl(k,ki)=odl(k,ki+iishift)
3872                sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
3873                l_homo(k,ki)=l_homo(k,ki+iishift)
3874               enddo
3875              enddo
3876              ii=ii-iishift
3877              lim_odl=lim_odl-iishift
3878           endif
3879          enddo
3880         enddo
3881       endif
3882
3883       return
3884    10 stop "Error in fragment file"
3885       end