2 implicit real*8 (a-h,o-z)
8 include 'COMMON.CONTROL'
9 include 'COMMON.SBRIDGE'
10 include 'COMMON.IOUNITS'
12 C Read force-field parameters except weights
14 C Read job setup parameters
16 C Read control parameters for energy minimzation if required
17 if (minim) call read_minim
18 C Read MCM control parameters if required
19 if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
20 C Read MD control parameters if reqjuired
21 if (modecalc.eq.12) call read_MDpar
22 C Read MREMD control parameters if required
23 if (modecalc.eq.14) then
27 C Read MUCA control parameters if required
28 if (lmuca) call read_muca
29 C Read CSA control parameters if required (from fort.40 if exists
30 C otherwise from general input file)
31 csa if (modecalc.eq.8) then
32 csa inquire (file="fort.40",exist=file_exist)
33 csa if (.not.file_exist) call csaread
35 cfmc if (modecalc.eq.10) call mcmfread
36 C Read molecule information, molecule geometry, energy-term weights, and
37 C restraints if requested
39 C Print restraint information
41 if (.not. out1file .or. me.eq.king) then
44 &write (iout,'(a,i5,a)') "The following",nhpb-nss,
45 & " distance constraints have been imposed"
47 write (iout,'(3i6,i2,3f10.5)') i-nss,ihpb(i),jhpb(i),
48 & ibecarb(i),dhpb(i),dhpb1(i),forcon(i)
53 c print *,"Processor",myrank," leaves READRTNS"
56 C-------------------------------------------------------------------------------
57 subroutine read_control
61 implicit real*8 (a-h,o-z)
65 logical OKRandom, prng_restart
68 include 'COMMON.IOUNITS'
69 include 'COMMON.TIME1'
70 include 'COMMON.THREAD'
71 include 'COMMON.SBRIDGE'
72 include 'COMMON.CONTROL'
75 include 'COMMON.HEADER'
76 csa include 'COMMON.CSA'
77 include 'COMMON.CHAIN'
80 include 'COMMON.FFIELD'
81 include 'COMMON.SETUP'
82 COMMON /MACHSW/ KDIAG,ICORFL,IXDR
83 character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
85 character*320 controlcard
90 read (INP,'(a)') titel
91 call card_concat(controlcard)
92 c out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
93 c print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
94 call reada(controlcard,'SEED',seed,0.0D0)
95 call random_init(seed)
96 C Set up the time limit (caution! The time must be input in minutes!)
97 read_cart=index(controlcard,'READ_CART').gt.0
98 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
99 call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
100 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
101 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
102 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
103 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
104 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
105 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
106 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
107 call reada(controlcard,'DRMS',drms,0.1D0)
108 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
109 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
110 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
111 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
112 write (iout,'(a,f10.1)')'DRMS = ',drms
113 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
114 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
116 call readi(controlcard,'NZ_START',nz_start,0)
117 call readi(controlcard,'NZ_END',nz_end,0)
118 call readi(controlcard,'IZ_SC',iz_sc,0)
120 safety = 60.0d0*safety
123 call reada(controlcard,"T_BATH",t_bath,300.0d0)
124 minim=(index(controlcard,'MINIMIZE').gt.0)
125 dccart=(index(controlcard,'CART').gt.0)
126 overlapsc=(index(controlcard,'OVERLAP').gt.0)
127 overlapsc=.not.overlapsc
128 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
129 searchsc=.not.searchsc
130 sideadd=(index(controlcard,'SIDEADD').gt.0)
131 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
132 outpdb=(index(controlcard,'PDBOUT').gt.0)
133 outmol2=(index(controlcard,'MOL2OUT').gt.0)
134 pdbref=(index(controlcard,'PDBREF').gt.0)
135 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
136 indpdb=index(controlcard,'PDBSTART')
137 extconf=(index(controlcard,'EXTCONF').gt.0)
138 call readi(controlcard,'IPRINT',iprint,0)
139 call readi(controlcard,'MAXGEN',maxgen,10000)
140 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
141 call readi(controlcard,"KDIAG",kdiag,0)
142 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
143 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
144 & write (iout,*) "RESCALE_MODE",rescale_mode
145 split_ene=index(controlcard,'SPLIT_ENE').gt.0
146 if (index(controlcard,'REGULAR').gt.0.0D0) then
147 call reada(controlcard,'WEIDIS',weidis,0.1D0)
151 if (index(controlcard,'CHECKGRAD').gt.0) then
153 if (index(controlcard,'CART').gt.0) then
155 elseif (index(controlcard,'CARINT').gt.0) then
160 elseif (index(controlcard,'THREAD').gt.0) then
162 call readi(controlcard,'THREAD',nthread,0)
163 if (nthread.gt.0) then
164 call reada(controlcard,'WEIDIS',weidis,0.1D0)
167 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
168 stop 'Error termination in Read_Control.'
170 else if (index(controlcard,'MCMA').gt.0) then
172 else if (index(controlcard,'MCEE').gt.0) then
174 else if (index(controlcard,'MULTCONF').gt.0) then
176 else if (index(controlcard,'MAP').gt.0) then
178 call readi(controlcard,'MAP',nmap,0)
179 else if (index(controlcard,'CSA').gt.0) then
180 write(*,*) "CSA not supported in this version"
183 crc else if (index(controlcard,'ZSCORE').gt.0) then
185 crc ZSCORE is rm from UNRES, modecalc=9 is available
188 cfcm else if (index(controlcard,'MCMF').gt.0) then
190 else if (index(controlcard,'SOFTREG').gt.0) then
192 else if (index(controlcard,'CHECK_BOND').gt.0) then
194 else if (index(controlcard,'TEST').gt.0) then
196 else if (index(controlcard,'MD').gt.0) then
198 else if (index(controlcard,'RE ').gt.0) then
202 lmuca=index(controlcard,'MUCA').gt.0
203 call readi(controlcard,'MUCADYN',mucadyn,0)
204 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
205 if (lmuca .and. (me.eq.king .or. .not.out1file ))
207 write (iout,*) 'MUCADYN=',mucadyn
208 write (iout,*) 'MUCASMOOTH=',muca_smooth
211 iscode=index(controlcard,'ONE_LETTER')
212 indphi=index(controlcard,'PHI')
213 indback=index(controlcard,'BACK')
214 iranconf=index(controlcard,'RAND_CONF')
215 i2ndstr=index(controlcard,'USE_SEC_PRED')
216 gradout=index(controlcard,'GRADOUT').gt.0
217 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
219 if(me.eq.king.or..not.out1file)
220 & write (iout,'(2a)') diagmeth(kdiag),
221 & ' routine used to diagonalize matrices.'
224 c--------------------------------------------------------------------------
225 subroutine read_REMDpar
229 implicit real*8 (a-h,o-z)
231 include 'COMMON.IOUNITS'
232 include 'COMMON.TIME1'
235 include 'COMMON.LANGEVIN'
237 include 'COMMON.LANGEVIN.lang0'
239 include 'COMMON.INTERACT'
240 include 'COMMON.NAMES'
242 include 'COMMON.REMD'
243 include 'COMMON.CONTROL'
244 include 'COMMON.SETUP'
246 character*320 controlcard
247 character*3200 controlcard1
248 integer iremd_m_total
250 if(me.eq.king.or..not.out1file)
251 & write (iout,*) "REMD setup"
253 call card_concat(controlcard)
254 call readi(controlcard,"NREP",nrep,3)
255 call readi(controlcard,"NSTEX",nstex,1000)
256 call reada(controlcard,"RETMIN",retmin,10.0d0)
257 call reada(controlcard,"RETMAX",retmax,1000.0d0)
258 mremdsync=(index(controlcard,'SYNC').gt.0)
259 call readi(controlcard,"NSYN",i_sync_step,100)
260 restart1file=(index(controlcard,'REST1FILE').gt.0)
261 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
262 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
263 if(max_cache_traj_use.gt.max_cache_traj)
264 & max_cache_traj_use=max_cache_traj
265 if(me.eq.king.or..not.out1file) then
266 cd if (traj1file) then
267 crc caching is in testing - NTWX is not ignored
268 cd write (iout,*) "NTWX value is ignored"
269 cd write (iout,*) " trajectory is stored to one file by master"
270 cd write (iout,*) " before exchange at NSTEX intervals"
272 write (iout,*) "NREP= ",nrep
273 write (iout,*) "NSTEX= ",nstex
274 write (iout,*) "SYNC= ",mremdsync
275 write (iout,*) "NSYN= ",i_sync_step
276 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
279 t_exchange_only=(index(controlcard,'TONLY').gt.0)
280 call readi(controlcard,"HREMD",hremd,0)
281 if((me.eq.king.or..not.out1file).and.hremd.gt.0) then
282 write (iout,*) "Hamiltonian REMD with ",hremd," sets of weights"
284 if(usampl.and.hremd.gt.0) then
286 & "========== ERROR: USAMPL and HREMD cannot be used together"
288 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
295 if (index(controlcard,'TLIST').gt.0) then
297 call card_concat(controlcard1)
298 read(controlcard1,*) (remd_t(i),i=1,nrep)
299 if(me.eq.king.or..not.out1file)
300 & write (iout,*)'tlist',(remd_t(i),i=1,nrep)
303 if (index(controlcard,'MLIST').gt.0) then
305 call card_concat(controlcard1)
306 read(controlcard1,*) (remd_m(i),i=1,nrep)
307 if(me.eq.king.or..not.out1file) then
308 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
311 iremd_m_total=iremd_m_total+remd_m(i)
314 write (iout,*) 'Total number of replicas ',
315 & iremd_m_total*hremd
317 write (iout,*) 'Total number of replicas ',iremd_m_total
321 if(me.eq.king.or..not.out1file)
322 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
325 c--------------------------------------------------------------------------
326 subroutine read_MDpar
330 implicit real*8 (a-h,o-z)
332 include 'COMMON.IOUNITS'
333 include 'COMMON.TIME1'
336 include 'COMMON.LANGEVIN'
338 include 'COMMON.LANGEVIN.lang0'
340 include 'COMMON.INTERACT'
341 include 'COMMON.NAMES'
343 include 'COMMON.SETUP'
344 include 'COMMON.CONTROL'
345 include 'COMMON.SPLITELE'
347 character*320 controlcard
349 call card_concat(controlcard)
350 call readi(controlcard,"NSTEP",n_timestep,1000000)
351 call readi(controlcard,"NTWE",ntwe,100)
352 call readi(controlcard,"NTWX",ntwx,1000)
353 call reada(controlcard,"DT",d_time,1.0d-1)
354 call reada(controlcard,"DVMAX",dvmax,2.0d1)
355 call reada(controlcard,"DAMAX",damax,1.0d1)
356 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
357 call readi(controlcard,"LANG",lang,0)
358 RESPA = index(controlcard,"RESPA") .gt. 0
359 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
360 ntime_split0=ntime_split
361 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
362 ntime_split0=ntime_split
363 call reada(controlcard,"R_CUT",r_cut,2.0d0)
364 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
365 rest = index(controlcard,"REST").gt.0
366 tbf = index(controlcard,"TBF").gt.0
367 call readi(controlcard,"HMC",hmc,0)
368 tnp = index(controlcard,"NOSEPOINCARE99").gt.0
369 tnp1 = index(controlcard,"NOSEPOINCARE01").gt.0
370 tnh = index(controlcard,"NOSEHOOVER96").gt.0
371 if (RESPA.and.tnh)then
372 xiresp = index(controlcard,"XIRESP").gt.0
374 call reada(controlcard,"Q_NP",Q_np,0.1d0)
375 usampl = index(controlcard,"USAMPL").gt.0
377 mdpdb = index(controlcard,"MDPDB").gt.0
378 call reada(controlcard,"T_BATH",t_bath,300.0d0)
379 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
380 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
381 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
382 if (count_reset_moment.eq.0) count_reset_moment=1000000000
383 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
384 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
385 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
386 if (count_reset_vel.eq.0) count_reset_vel=1000000000
387 large = index(controlcard,"LARGE").gt.0
388 print_compon = index(controlcard,"PRINT_COMPON").gt.0
389 rattle = index(controlcard,"RATTLE").gt.0
390 c if performing umbrella sampling, fragments constrained are read from the fragment file
396 if(me.eq.king.or..not.out1file) then
398 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
400 write (iout,'(a)') "The units are:"
401 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
402 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
403 & " acceleration: angstrom/(48.9 fs)**2"
404 write (iout,'(a)') "energy: kcal/mol, temperature: K"
406 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
407 write (iout,'(a60,f10.5,a)')
408 & "Initial time step of numerical integration:",d_time,
410 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
412 write (iout,'(2a,i4,a)')
413 & "A-MTS algorithm used; initial time step for fast-varying",
414 & " short-range forces split into",ntime_split," steps."
415 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
416 & r_cut," lambda",rlamb
418 write (iout,'(2a,f10.5)')
419 & "Maximum acceleration threshold to reduce the time step",
420 & "/increase split number:",damax
421 write (iout,'(2a,f10.5)')
422 & "Maximum predicted energy drift to reduce the timestep",
423 & "/increase split number:",edriftmax
424 write (iout,'(a60,f10.5)')
425 & "Maximum velocity threshold to reduce velocities:",dvmax
426 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
427 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
428 if (rattle) write (iout,'(a60)')
429 & "Rattle algorithm used to constrain the virtual bonds"
433 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
434 call reada(controlcard,"RWAT",rwat,1.4d0)
435 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
436 surfarea=index(controlcard,"SURFAREA").gt.0
437 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
438 if(me.eq.king.or..not.out1file)then
439 write (iout,'(/a,$)') "Langevin dynamics calculation"
442 & " with direct integration of Langevin equations"
443 else if (lang.eq.2) then
444 write (iout,'(a/)') " with TINKER stochasic MD integrator"
445 else if (lang.eq.3) then
446 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
447 else if (lang.eq.4) then
448 write (iout,'(a/)') " in overdamped mode"
450 write (iout,'(//a,i5)')
451 & "=========== ERROR: Unknown Langevin dynamics mode:",lang
454 write (iout,'(a60,f10.5)') "Temperature:",t_bath
455 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
456 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
457 write (iout,'(a60,f10.5)')
458 & "Scaling factor of the friction forces:",scal_fric
459 if (surfarea) write (iout,'(2a,i10,a)')
460 & "Friction coefficients will be scaled by solvent-accessible",
461 & " surface area every",reset_fricmat," steps."
463 c Calculate friction coefficients and bounds of stochastic forces
464 eta=6*pi*cPoise*etawat
465 if(me.eq.king.or..not.out1file)
466 & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
468 gamp=scal_fric*(pstok+rwat)*eta
469 stdfp=dsqrt(2*Rb*t_bath/d_time)
471 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
472 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
474 if(me.eq.king.or..not.out1file)then
475 write (iout,'(/2a/)')
476 & "Radii of site types and friction coefficients and std's of",
477 & " stochastic forces of fully exposed sites"
478 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
480 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
481 & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
485 if(me.eq.king.or..not.out1file)then
486 write (iout,'(a)') "Berendsen bath calculation"
487 write (iout,'(a60,f10.5)') "Temperature:",t_bath
488 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
490 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
491 & count_reset_moment," steps"
493 & write (iout,'(a,i10,a)')
494 & "Velocities will be reset at random every",count_reset_vel,
497 else if (tnp .or. tnp1 .or. tnh) then
498 if (tnp .or. tnp1) then
499 write (iout,'(a)') "Nose-Poincare bath calculation"
500 if (tnp) write (iout,'(a)')
501 & "J.Comput.Phys. 151 114 (1999) S.D.Bond B.J.Leimkuhler B.B.Laird"
502 if (tnp1) write (iout,'(a)') "JPSJ 70 75 (2001) S. Nose"
504 write (iout,'(a)') "Nose-Hoover bath calculation"
505 write (iout,'(a)') "Mol.Phys. 87 1117 (1996) Martyna et al."
515 WDTI(i) = 1.0*d_time/nresn
522 write (iout,'(a)') "NVT-XI-RESPA algorithm"
524 write (iout,'(a)') "NVT-XO-RESPA algorithm"
527 WDTIi(i) = 1.0*d_time/nresn/ntime_split
535 write (iout,'(a60,f10.5)') "Temperature:",t_bath
536 write (iout,'(a60,f10.5)') "Q =",Q_np
538 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
539 & count_reset_moment," steps"
541 & write (iout,'(a,i10,a)')
542 & "Velocities will be reset at random every",count_reset_vel,
545 else if (hmc.gt.0) then
546 write (iout,'(a)') "Hybrid Monte Carlo calculation"
547 write (iout,'(a60,f10.5)') "Temperature:",t_bath
548 write (iout,'(a60,i10)')
549 & "Number of MD steps between Metropolis tests:",hmc
552 if(me.eq.king.or..not.out1file)
553 & write (iout,'(a31)') "Microcanonical mode calculation"
555 if(me.eq.king.or..not.out1file)then
556 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
558 write(iout,*) "MD running with constraints."
559 write(iout,*) "Equilibration time ", eq_time, " mtus."
560 write(iout,*) "Constraining ", nfrag," fragments."
561 write(iout,*) "Length of each fragment, weight and q0:"
563 write (iout,*) "Set of restraints #",iset
565 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
566 & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
568 write(iout,*) "constraints between ", npair, "fragments."
569 write(iout,*) "constraint pairs, weights and q0:"
571 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
572 & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
574 write(iout,*) "angle constraints within ", nfrag_back,
575 & "backbone fragments."
576 write(iout,*) "fragment, weights:"
578 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
579 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
580 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
583 iset=mod(kolor,nset)+1
586 if(me.eq.king.or..not.out1file)
587 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
590 c------------------------------------------------------------------------------
593 C Read molecular data.
595 implicit real*8 (a-h,o-z)
601 include 'COMMON.IOUNITS'
604 include 'COMMON.INTERACT'
605 include 'COMMON.LOCAL'
606 include 'COMMON.NAMES'
607 include 'COMMON.CHAIN'
608 include 'COMMON.FFIELD'
609 include 'COMMON.SBRIDGE'
610 include 'COMMON.HEADER'
611 include 'COMMON.CONTROL'
612 include 'COMMON.DBASE'
613 include 'COMMON.THREAD'
614 include 'COMMON.CONTACTS'
615 include 'COMMON.TORCNSTR'
616 include 'COMMON.TIME1'
617 include 'COMMON.BOUNDS'
619 include 'COMMON.REMD'
620 include 'COMMON.SETUP'
621 character*4 sequence(maxres)
623 double precision x(maxvar)
624 character*256 pdbfile
625 character*320 weightcard
626 character*80 weightcard_t,ucase
627 dimension itype_pdb(maxres)
628 common /pizda/ itype_pdb
629 logical seq_comp,fail
630 double precision energia(0:n_ene)
636 C Read weights of the subsequent energy terms.
649 if(me.eq.king.or..not.out1file) then
650 write (iout,*) 'Reading ',hremd,' sets of weights for HREMD'
651 write (iout,*) 'Current weights for processor ',
652 & me,' set ',i2set(me)
656 call card_concat(weightcard)
657 call reada(weightcard,'WLONG',wlong,1.0D0)
658 call reada(weightcard,'WSC',wsc,wlong)
659 call reada(weightcard,'WSCP',wscp,wlong)
660 call reada(weightcard,'WELEC',welec,1.0D0)
661 call reada(weightcard,'WVDWPP',wvdwpp,welec)
662 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
663 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
664 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
665 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
666 call reada(weightcard,'WTURN3',wturn3,1.0D0)
667 call reada(weightcard,'WTURN4',wturn4,1.0D0)
668 call reada(weightcard,'WTURN6',wturn6,1.0D0)
669 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
670 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
671 call reada(weightcard,'WBOND',wbond,1.0D0)
672 call reada(weightcard,'WTOR',wtor,1.0D0)
673 call reada(weightcard,'WTORD',wtor_d,1.0D0)
674 call reada(weightcard,'WANG',wang,1.0D0)
675 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
676 call reada(weightcard,'SCAL14',scal14,0.4D0)
677 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
678 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
679 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
680 call reada(weightcard,'TEMP0',temp0,300.0d0)
681 if (index(weightcard,'SOFT').gt.0) ipot=6
682 C 12/1/95 Added weight for the multi-body term WCORR
683 call reada(weightcard,'WCORRH',wcorr,1.0D0)
684 if (wcorr4.gt.0.0d0) wcorr=wcorr4
692 hweights(i,7)=wel_loc
695 hweights(i,10)=wturn6
697 hweights(i,12)=wscloc
699 hweights(i,14)=wtor_d
700 hweights(i,15)=wstrain
701 hweights(i,16)=wvdwpp
703 hweights(i,18)=scal14
704 hweights(i,21)=wsccor
709 weights(i)=hweights(i2set(me),i)
733 call card_concat(weightcard)
734 call reada(weightcard,'WLONG',wlong,1.0D0)
735 call reada(weightcard,'WSC',wsc,wlong)
736 call reada(weightcard,'WSCP',wscp,wlong)
737 call reada(weightcard,'WELEC',welec,1.0D0)
738 call reada(weightcard,'WVDWPP',wvdwpp,welec)
739 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
740 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
741 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
742 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
743 call reada(weightcard,'WTURN3',wturn3,1.0D0)
744 call reada(weightcard,'WTURN4',wturn4,1.0D0)
745 call reada(weightcard,'WTURN6',wturn6,1.0D0)
746 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
747 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
748 call reada(weightcard,'WBOND',wbond,1.0D0)
749 call reada(weightcard,'WTOR',wtor,1.0D0)
750 call reada(weightcard,'WTORD',wtor_d,1.0D0)
751 call reada(weightcard,'WANG',wang,1.0D0)
752 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
753 call reada(weightcard,'SCAL14',scal14,0.4D0)
754 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
755 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
756 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
757 call reada(weightcard,'TEMP0',temp0,300.0d0)
758 if (index(weightcard,'SOFT').gt.0) ipot=6
759 C 12/1/95 Added weight for the multi-body term WCORR
760 call reada(weightcard,'WCORRH',wcorr,1.0D0)
761 if (wcorr4.gt.0.0d0) wcorr=wcorr4
783 if(me.eq.king.or..not.out1file)
784 & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
785 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
787 10 format (/'Energy-term weights (unscaled):'//
788 & 'WSCC= ',f10.6,' (SC-SC)'/
789 & 'WSCP= ',f10.6,' (SC-p)'/
790 & 'WELEC= ',f10.6,' (p-p electr)'/
791 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
792 & 'WBOND= ',f10.6,' (stretching)'/
793 & 'WANG= ',f10.6,' (bending)'/
794 & 'WSCLOC= ',f10.6,' (SC local)'/
795 & 'WTOR= ',f10.6,' (torsional)'/
796 & 'WTORD= ',f10.6,' (double torsional)'/
797 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
798 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
799 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
800 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
801 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
802 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
803 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
804 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
805 & 'WTURN6= ',f10.6,' (turns, 6th order)')
806 if(me.eq.king.or..not.out1file)then
807 if (wcorr4.gt.0.0d0) then
808 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
809 & 'between contact pairs of peptide groups'
810 write (iout,'(2(a,f5.3/))')
811 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
812 & 'Range of quenching the correlation terms:',2*delt_corr
813 else if (wcorr.gt.0.0d0) then
814 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
815 & 'between contact pairs of peptide groups'
817 write (iout,'(a,f8.3)')
818 & 'Scaling factor of 1,4 SC-p interactions:',scal14
819 write (iout,'(a,f8.3)')
820 & 'General scaling factor of SC-p interactions:',scalscp
822 r0_corr=cutoff_corr-delt_corr
824 aad(i,1)=scalscp*aad(i,1)
825 aad(i,2)=scalscp*aad(i,2)
826 bad(i,1)=scalscp*bad(i,1)
827 bad(i,2)=scalscp*bad(i,2)
829 call rescale_weights(t_bath)
830 if(me.eq.king.or..not.out1file)
831 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
832 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
834 22 format (/'Energy-term weights (scaled):'//
835 & 'WSCC= ',f10.6,' (SC-SC)'/
836 & 'WSCP= ',f10.6,' (SC-p)'/
837 & 'WELEC= ',f10.6,' (p-p electr)'/
838 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
839 & 'WBOND= ',f10.6,' (stretching)'/
840 & 'WANG= ',f10.6,' (bending)'/
841 & 'WSCLOC= ',f10.6,' (SC local)'/
842 & 'WTOR= ',f10.6,' (torsional)'/
843 & 'WTORD= ',f10.6,' (double torsional)'/
844 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
845 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
846 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
847 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
848 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
849 & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
850 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
851 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
852 & 'WTURN6= ',f10.6,' (turns, 6th order)')
853 if(me.eq.king.or..not.out1file)
854 & write (iout,*) "Reference temperature for weights calculation:",
856 call reada(weightcard,"D0CM",d0cm,3.78d0)
857 call reada(weightcard,"AKCM",akcm,15.1d0)
858 call reada(weightcard,"AKTH",akth,11.0d0)
859 call reada(weightcard,"AKCT",akct,12.0d0)
860 call reada(weightcard,"V1SS",v1ss,-1.08d0)
861 call reada(weightcard,"V2SS",v2ss,7.61d0)
862 call reada(weightcard,"V3SS",v3ss,13.7d0)
863 call reada(weightcard,"EBR",ebr,-5.50D0)
864 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
866 dyn_ss_mask(i)=.false.
870 dyn_ssbond_ij(i,j)=1.0d300
873 call reada(weightcard,"HT",Ht,0.0D0)
875 ss_depth=ebr/wsc-0.25*eps(1,1)
876 Ht=Ht/wsc-0.25*eps(1,1)
877 akcm=akcm*wstrain/wsc
878 akth=akth*wstrain/wsc
879 akct=akct*wstrain/wsc
880 v1ss=v1ss*wstrain/wsc
881 v2ss=v2ss*wstrain/wsc
882 v3ss=v3ss*wstrain/wsc
884 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
887 if(me.eq.king.or..not.out1file) then
888 write (iout,*) "Parameters of the SS-bond potential:"
889 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
891 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
892 write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
893 write (iout,*)" HT",Ht
894 print *,'indpdb=',indpdb,' pdbref=',pdbref
896 if (indpdb.gt.0 .or. pdbref) then
897 read(inp,'(a)') pdbfile
898 if(me.eq.king.or..not.out1file)
899 & write (iout,'(2a)') 'PDB data will be read from file ',
900 & pdbfile(:ilen(pdbfile))
901 open(ipdbin,file=pdbfile,status='old',err=33)
903 33 write (iout,'(a)') 'Error opening PDB file.'
906 c print *,'Begin reading pdb data'
908 c print *,'Finished reading pdb data'
909 if(me.eq.king.or..not.out1file)
910 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
911 & ' nstart_sup=',nstart_sup
913 itype_pdb(i)=itype(i)
917 nct=nstart_sup+nsup-1
918 call contact(.false.,ncont_ref,icont_ref,co)
921 C Following 2 lines for diagnostics; comment out if not needed
922 write (iout,*) "Before sideadd"
924 if(me.eq.king.or..not.out1file)
925 & write(iout,*)'Adding sidechains'
932 do while (fail.and.nsi.le.maxsi)
933 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
936 if(fail) write(iout,*)'Adding sidechain failed for res ',
937 & i,' after ',nsi,' trials'
940 C 10/03/12 Adam: Recalculate coordinates with new side chain positions
943 C Following 2 lines for diagnostics; comment out if not needed
944 c write (iout,*) "After sideadd"
947 if (indpdb.eq.0) then
948 C Read sequence if not taken from the pdb file.
950 c print *,'nres=',nres
951 if (iscode.gt.0) then
952 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
954 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
956 C Convert sequence to numeric code
958 itype(i)=rescode(i,sequence(i),iscode)
960 C Assign initial virtual bond lengths
966 vbld(i+nres)=dsc(itype(i))
967 vbld_inv(i+nres)=dsc_inv(itype(i))
968 c write (iout,*) "i",i," itype",itype(i),
969 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
973 c print '(20i4)',(itype(i),i=1,nres)
976 if (itype(i).eq.21 .or. itype(i+1).eq.21) then
978 if (itype(i).eq.21) then
982 else if (itype(i+1).ne.20) then
984 else if (itype(i).ne.20) then
991 if(me.eq.king.or..not.out1file)then
992 write (iout,*) "ITEL"
994 write (iout,*) i,itype(i),itel(i)
996 print *,'Call Read_Bridge.'
999 C 8/13/98 Set limits to generating the dihedral angles
1004 read (inp,*) ndih_constr
1005 if (ndih_constr.gt.0) then
1007 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
1008 if(me.eq.king.or..not.out1file)then
1010 & 'There are',ndih_constr,' constraints on phi angles.'
1012 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
1016 phi0(i)=deg2rad*phi0(i)
1017 drange(i)=deg2rad*drange(i)
1019 if(me.eq.king.or..not.out1file)
1020 & write (iout,*) 'FTORS',ftors
1023 phibound(1,ii) = phi0(i)-drange(i)
1024 phibound(2,ii) = phi0(i)+drange(i)
1029 if (me.eq.king) then
1031 write (iout,'(a)') 'Boundaries in phi angle sampling:'
1033 write (iout,'(a3,i5,2f10.1)')
1034 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
1040 cd print *,'NNT=',NNT,' NCT=',NCT
1041 if (itype(1).eq.21) nnt=2
1042 if (itype(nres).eq.21) nct=nct-1
1044 if(me.eq.king.or..not.out1file)
1045 & write (iout,'(a,i3)') 'nsup=',nsup
1047 if (nsup.le.(nct-nnt+1)) then
1048 do i=0,nct-nnt+1-nsup
1049 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
1055 & 'Error - sequences to be superposed do not match.'
1058 do i=0,nsup-(nct-nnt+1)
1059 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
1061 nstart_sup=nstart_sup+i
1067 & 'Error - sequences to be superposed do not match.'
1070 if (nsup.eq.0) nsup=nct-nnt
1071 if (nstart_sup.eq.0) nstart_sup=nnt
1072 if (nstart_seq.eq.0) nstart_seq=nnt
1073 if(me.eq.king.or..not.out1file)
1074 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
1075 & ' nstart_seq=',nstart_seq
1077 c--- Zscore rms -------
1078 if (nz_start.eq.0) nz_start=nnt
1079 if (nz_end.eq.0 .and. nsup.gt.0) then
1081 else if (nz_end.eq.0) then
1084 if(me.eq.king.or..not.out1file)then
1085 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1086 write (iout,*) 'IZ_SC=',iz_sc
1088 c----------------------
1091 if (.not.pdbref) then
1092 call read_angles(inp,*38)
1094 38 write (iout,'(a)') 'Error reading reference structure.'
1096 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1097 stop 'Error reading reference structure'
1101 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
1110 call contact(.true.,ncont_ref,icont_ref,co)
1112 if(me.eq.king.or..not.out1file)
1113 & write (iout,*) 'Contact order:',co
1115 if(me.eq.king.or..not.out1file)
1116 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1119 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1121 if(me.eq.king.or..not.out1file)
1122 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1123 & icont_ref(1,i),' ',
1124 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1128 c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1129 if (constr_dist.gt.0) then
1130 call read_dist_constr
1134 if (constr_homology.gt.0) then
1135 call read_constr_homology
1139 if (nhpb.gt.0) call hpb_partition
1140 c write (iout,*) "After read_dist_constr nhpb",nhpb
1142 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1143 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1144 & modecalc.ne.10) then
1145 C If input structure hasn't been supplied from the PDB file read or generate
1147 if (iranconf.eq.0 .and. .not. extconf) then
1148 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1149 & write (iout,'(a)') 'Initial geometry will be read in.'
1151 read(inp,'(8f10.5)',end=36,err=36)
1152 & ((c(l,k),l=1,3),k=1,nres),
1153 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1154 call int_from_cart1(.false.)
1157 dc(j,i)=c(j,i+1)-c(j,i)
1158 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1162 if (itype(i).ne.10) then
1164 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1165 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1171 call read_angles(inp,*36)
1174 36 write (iout,'(a)') 'Error reading angle file.'
1176 call mpi_finalize( MPI_COMM_WORLD,IERR )
1178 stop 'Error reading angle file.'
1180 else if (extconf) then
1181 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1182 & write (iout,'(a)') 'Extended chain initial geometry.'
1184 theta(i)=90d0*deg2rad
1187 phi(i)=180d0*deg2rad
1190 alph(i)=110d0*deg2rad
1193 omeg(i)=-120d0*deg2rad
1196 if(me.eq.king.or..not.out1file)
1197 & write (iout,'(a)') 'Random-generated initial geometry.'
1201 if (me.eq.king .or. fg_rank.eq.0 .and. (
1202 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1206 call gen_rand_conf(itmp,*30)
1208 30 write (iout,*) 'Failed to generate random conformation',
1209 & ', itrial=',itrial
1210 write (*,*) 'Processor:',me,
1211 & ' Failed to generate random conformation',
1221 write (iout,'(a,i3,a)') 'Processor:',me,
1222 & ' error in generating random conformation.'
1223 write (*,'(a,i3,a)') 'Processor:',me,
1224 & ' error in generating random conformation.'
1227 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1234 elseif (modecalc.eq.4) then
1235 read (inp,'(a)') intinname
1236 open (intin,file=intinname,status='old',err=333)
1237 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1238 & write (iout,'(a)') 'intinname',intinname
1239 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1241 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1243 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1245 stop 'Error opening angle file.'
1249 C Generate distance constraints, if the PDB structure is to be regularized.
1250 if (nthread.gt.0) then
1251 call read_threadbase
1254 if (me.eq.king .or. .not. out1file)
1256 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1257 write (iout,'(/a,i3,a)')
1258 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1259 write (iout,'(20i4)') (iss(i),i=1,ns)
1261 write(iout,*)"Running with dynamic disulfide-bond formation"
1263 write (iout,'(/a/)') 'Pre-formed links are:'
1269 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1270 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1276 if (ns.gt.0.and.dyn_ss) then
1280 forcon(i-nss)=forcon(i)
1287 dyn_ss_mask(iss(i))=.true.
1290 if (i2ndstr.gt.0) call secstrp2dihc
1291 c call geom_to_var(nvar,x)
1292 c call etotal(energia(0))
1293 c call enerprint(energia(0))
1294 c call briefout(0,etot)
1296 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1297 cd write (iout,'(a)') 'Variable list:'
1298 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1300 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1301 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1302 & 'Processor',myrank,': end reading molecular data.'
1306 c--------------------------------------------------------------------------
1307 logical function seq_comp(itypea,itypeb,length)
1309 integer length,itypea(length),itypeb(length)
1312 if (itypea(i).ne.itypeb(i)) then
1320 c-----------------------------------------------------------------------------
1321 subroutine read_bridge
1322 C Read information about disulfide bridges.
1323 implicit real*8 (a-h,o-z)
1324 include 'DIMENSIONS'
1328 include 'COMMON.IOUNITS'
1329 include 'COMMON.GEO'
1330 include 'COMMON.VAR'
1331 include 'COMMON.INTERACT'
1332 include 'COMMON.LOCAL'
1333 include 'COMMON.NAMES'
1334 include 'COMMON.CHAIN'
1335 include 'COMMON.FFIELD'
1336 include 'COMMON.SBRIDGE'
1337 include 'COMMON.HEADER'
1338 include 'COMMON.CONTROL'
1339 include 'COMMON.DBASE'
1340 include 'COMMON.THREAD'
1341 include 'COMMON.TIME1'
1342 include 'COMMON.SETUP'
1343 C Read bridging residues.
1344 read (inp,*) ns,(iss(i),i=1,ns)
1346 if(me.eq.king.or..not.out1file)
1347 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1348 C Check whether the specified bridging residues are cystines.
1350 if (itype(iss(i)).ne.1) then
1351 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1352 & 'Do you REALLY think that the residue ',
1353 & restyp(itype(iss(i))),i,
1354 & ' can form a disulfide bridge?!!!'
1355 write (*,'(2a,i3,a)')
1356 & 'Do you REALLY think that the residue ',
1357 & restyp(itype(iss(i))),i,
1358 & ' can form a disulfide bridge?!!!'
1360 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1365 C Read preformed bridges.
1367 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1369 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1372 C Check if the residues involved in bridges are in the specified list of
1373 C bridging residues.
1376 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1377 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1378 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1379 & ' contains residues present in other pairs.'
1380 write (*,'(a,i3,a)') 'Disulfide pair',i,
1381 & ' contains residues present in other pairs.'
1383 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1389 if (ihpb(i).eq.iss(j)) goto 10
1391 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1394 if (jhpb(i).eq.iss(j)) goto 20
1396 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1402 ihpb(i)=ihpb(i)+nres
1403 jhpb(i)=jhpb(i)+nres
1409 c----------------------------------------------------------------------------
1410 subroutine read_x(kanal,*)
1411 implicit real*8 (a-h,o-z)
1412 include 'DIMENSIONS'
1413 include 'COMMON.GEO'
1414 include 'COMMON.VAR'
1415 include 'COMMON.CHAIN'
1416 include 'COMMON.IOUNITS'
1417 include 'COMMON.CONTROL'
1418 include 'COMMON.LOCAL'
1419 include 'COMMON.INTERACT'
1420 c Read coordinates from input
1422 read(kanal,'(8f10.5)',end=10,err=10)
1423 & ((c(l,k),l=1,3),k=1,nres),
1424 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1427 c(j,2*nres)=c(j,nres)
1429 call int_from_cart1(.false.)
1432 dc(j,i)=c(j,i+1)-c(j,i)
1433 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1437 if (itype(i).ne.10) then
1439 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1440 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1448 c----------------------------------------------------------------------------
1449 subroutine read_threadbase
1450 implicit real*8 (a-h,o-z)
1451 include 'DIMENSIONS'
1452 include 'COMMON.IOUNITS'
1453 include 'COMMON.GEO'
1454 include 'COMMON.VAR'
1455 include 'COMMON.INTERACT'
1456 include 'COMMON.LOCAL'
1457 include 'COMMON.NAMES'
1458 include 'COMMON.CHAIN'
1459 include 'COMMON.FFIELD'
1460 include 'COMMON.SBRIDGE'
1461 include 'COMMON.HEADER'
1462 include 'COMMON.CONTROL'
1463 include 'COMMON.DBASE'
1464 include 'COMMON.THREAD'
1465 include 'COMMON.TIME1'
1466 C Read pattern database for threading.
1467 read (icbase,*) nseq
1469 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1470 & nres_base(2,i),nres_base(3,i)
1471 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1473 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1474 c & nres_base(2,i),nres_base(3,i)
1475 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1479 if (weidis.eq.0.0D0) weidis=0.1D0
1488 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1489 write (iout,'(a,i5)') 'nexcl: ',nexcl
1490 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1493 c------------------------------------------------------------------------------
1494 subroutine setup_var
1495 implicit real*8 (a-h,o-z)
1496 include 'DIMENSIONS'
1497 include 'COMMON.IOUNITS'
1498 include 'COMMON.GEO'
1499 include 'COMMON.VAR'
1500 include 'COMMON.INTERACT'
1501 include 'COMMON.LOCAL'
1502 include 'COMMON.NAMES'
1503 include 'COMMON.CHAIN'
1504 include 'COMMON.FFIELD'
1505 include 'COMMON.SBRIDGE'
1506 include 'COMMON.HEADER'
1507 include 'COMMON.CONTROL'
1508 include 'COMMON.DBASE'
1509 include 'COMMON.THREAD'
1510 include 'COMMON.TIME1'
1511 C Set up variable list.
1517 if (itype(i).ne.10) then
1519 ialph(i,1)=nvar+nside
1523 if (indphi.gt.0) then
1525 else if (indback.gt.0) then
1530 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1533 c----------------------------------------------------------------------------
1534 subroutine gen_dist_constr
1535 C Generate CA distance constraints.
1536 implicit real*8 (a-h,o-z)
1537 include 'DIMENSIONS'
1538 include 'COMMON.IOUNITS'
1539 include 'COMMON.GEO'
1540 include 'COMMON.VAR'
1541 include 'COMMON.INTERACT'
1542 include 'COMMON.LOCAL'
1543 include 'COMMON.NAMES'
1544 include 'COMMON.CHAIN'
1545 include 'COMMON.FFIELD'
1546 include 'COMMON.SBRIDGE'
1547 include 'COMMON.HEADER'
1548 include 'COMMON.CONTROL'
1549 include 'COMMON.DBASE'
1550 include 'COMMON.THREAD'
1551 include 'COMMON.TIME1'
1552 dimension itype_pdb(maxres)
1553 common /pizda/ itype_pdb
1555 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1556 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1557 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1559 do i=nstart_sup,nstart_sup+nsup-1
1560 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1561 cd & ' seq_pdb', restyp(itype_pdb(i))
1562 do j=i+2,nstart_sup+nsup-1
1564 ihpb(nhpb)=i+nstart_seq-nstart_sup
1565 jhpb(nhpb)=j+nstart_seq-nstart_sup
1567 dhpb(nhpb)=dist(i,j)
1570 cd write (iout,'(a)') 'Distance constraints:'
1575 cd if (ii.gt.nres) then
1580 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1581 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1582 cd & dhpb(i),forcon(i)
1586 c----------------------------------------------------------------------------
1588 implicit real*8 (a-h,o-z)
1589 include 'DIMENSIONS'
1590 include 'COMMON.MAP'
1591 include 'COMMON.IOUNITS'
1592 character*3 angid(4) /'THE','PHI','ALP','OME'/
1593 character*80 mapcard,ucase
1595 read (inp,'(a)') mapcard
1596 mapcard=ucase(mapcard)
1597 if (index(mapcard,'PHI').gt.0) then
1599 else if (index(mapcard,'THE').gt.0) then
1601 else if (index(mapcard,'ALP').gt.0) then
1603 else if (index(mapcard,'OME').gt.0) then
1606 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1607 stop 'Error - illegal variable spec in MAP card.'
1609 call readi (mapcard,'RES1',res1(imap),0)
1610 call readi (mapcard,'RES2',res2(imap),0)
1611 if (res1(imap).eq.0) then
1612 res1(imap)=res2(imap)
1613 else if (res2(imap).eq.0) then
1614 res2(imap)=res1(imap)
1616 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1618 & 'Error - illegal definition of variable group in MAP.'
1619 stop 'Error - illegal definition of variable group in MAP.'
1621 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1622 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1623 call readi(mapcard,'NSTEP',nstep(imap),0)
1624 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1626 & 'Illegal boundary and/or step size specification in MAP.'
1627 stop 'Illegal boundary and/or step size specification in MAP.'
1632 c----------------------------------------------------------------------------
1633 csa subroutine csaread
1634 csa implicit real*8 (a-h,o-z)
1635 csa include 'DIMENSIONS'
1636 csa include 'COMMON.IOUNITS'
1637 csa include 'COMMON.GEO'
1638 csa include 'COMMON.CSA'
1639 csa include 'COMMON.BANK'
1640 csa include 'COMMON.CONTROL'
1641 csa character*80 ucase
1642 csa character*620 mcmcard
1643 csa call card_concat(mcmcard)
1645 csa call readi(mcmcard,'NCONF',nconf,50)
1646 csa call readi(mcmcard,'NADD',nadd,0)
1647 csa call readi(mcmcard,'JSTART',jstart,1)
1648 csa call readi(mcmcard,'JEND',jend,1)
1649 csa call readi(mcmcard,'NSTMAX',nstmax,500000)
1650 csa call readi(mcmcard,'N0',n0,1)
1651 csa call readi(mcmcard,'N1',n1,6)
1652 csa call readi(mcmcard,'N2',n2,4)
1653 csa call readi(mcmcard,'N3',n3,0)
1654 csa call readi(mcmcard,'N4',n4,0)
1655 csa call readi(mcmcard,'N5',n5,0)
1656 csa call readi(mcmcard,'N6',n6,10)
1657 csa call readi(mcmcard,'N7',n7,0)
1658 csa call readi(mcmcard,'N8',n8,0)
1659 csa call readi(mcmcard,'N9',n9,0)
1660 csa call readi(mcmcard,'N14',n14,0)
1661 csa call readi(mcmcard,'N15',n15,0)
1662 csa call readi(mcmcard,'N16',n16,0)
1663 csa call readi(mcmcard,'N17',n17,0)
1664 csa call readi(mcmcard,'N18',n18,0)
1666 csa vdisulf=(index(mcmcard,'DYNSS').gt.0)
1668 csa call readi(mcmcard,'NDIFF',ndiff,2)
1669 csa call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1670 csa call readi(mcmcard,'IS1',is1,1)
1671 csa call readi(mcmcard,'IS2',is2,8)
1672 csa call readi(mcmcard,'NRAN0',nran0,4)
1673 csa call readi(mcmcard,'NRAN1',nran1,2)
1674 csa call readi(mcmcard,'IRR',irr,1)
1675 csa call readi(mcmcard,'NSEED',nseed,20)
1676 csa call readi(mcmcard,'NTOTAL',ntotal,10000)
1677 csa call reada(mcmcard,'CUT1',cut1,2.0d0)
1678 csa call reada(mcmcard,'CUT2',cut2,5.0d0)
1679 csa call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1680 csa call readi(mcmcard,'ICMAX',icmax,3)
1681 csa call readi(mcmcard,'IRESTART',irestart,0)
1682 csac!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1685 csa call reada(mcmcard,'DELE',dele,20.0d0)
1686 csa call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1687 csa call readi(mcmcard,'IREF',iref,0)
1688 csa call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1689 csa call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1690 csa call readi(mcmcard,'NCONF_IN',nconf_in,0)
1691 csa call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1692 csa write (iout,*) "NCONF_IN",nconf_in
1695 c----------------------------------------------------------------------------
1696 cfmc subroutine mcmfread
1697 cfmc implicit real*8 (a-h,o-z)
1698 cfmc include 'DIMENSIONS'
1699 cfmc include 'COMMON.MCMF'
1700 cfmc include 'COMMON.IOUNITS'
1701 cfmc include 'COMMON.GEO'
1702 cfmc character*80 ucase
1703 cfmc character*620 mcmcard
1704 cfmc call card_concat(mcmcard)
1706 cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
1707 cfmc write(iout,*)'MAXRANT=',maxrant
1708 cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1709 cfmc write(iout,*)'MAXFAM=',maxfam
1710 cfmc call readi(mcmcard,'NNET1',nnet1,5)
1711 cfmc write(iout,*)'NNET1=',nnet1
1712 cfmc call readi(mcmcard,'NNET2',nnet2,4)
1713 cfmc write(iout,*)'NNET2=',nnet2
1714 cfmc call readi(mcmcard,'NNET3',nnet3,4)
1715 cfmc write(iout,*)'NNET3=',nnet3
1716 cfmc call readi(mcmcard,'ILASTT',ilastt,0)
1717 cfmc write(iout,*)'ILASTT=',ilastt
1718 cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1719 cfmc write(iout,*)'MAXSTR=',maxstr
1720 cfmc maxstr_f=maxstr/maxfam
1721 cfmc write(iout,*)'MAXSTR_F=',maxstr_f
1722 cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
1723 cfmc write(iout,*)'NMCMF=',nmcmf
1724 cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1725 cfmc write(iout,*)'IFOCUS=',ifocus
1726 cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1727 cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
1728 cfmc call readi(mcmcard,'INTPRT',intprt,1000)
1729 cfmc write(iout,*)'INTPRT=',intprt
1730 cfmc call readi(mcmcard,'IPRT',iprt,100)
1731 cfmc write(iout,*)'IPRT=',iprt
1732 cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
1733 cfmc write(iout,*)'IMAXTR=',imaxtr
1734 cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
1735 cfmc write(iout,*)'MAXEVEN=',maxeven
1736 cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1737 cfmc write(iout,*)'MAXEVEN1=',maxeven1
1738 cfmc call readi(mcmcard,'INIMIN',inimin,200)
1739 cfmc write(iout,*)'INIMIN=',inimin
1740 cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1741 cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
1742 cfmc call readi(mcmcard,'NTHREAD',nthread,5)
1743 cfmc write(iout,*)'NTHREAD=',nthread
1744 cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1745 cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1746 cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
1747 cfmc write(iout,*)'MAXPERT=',maxpert
1748 cfmc call readi(mcmcard,'IRMSD',irmsd,1)
1749 cfmc write(iout,*)'IRMSD=',irmsd
1750 cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1751 cfmc write(iout,*)'DENEMIN=',denemin
1752 cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1753 cfmc write(iout,*)'RCUT1S=',rcut1s
1754 cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1755 cfmc write(iout,*)'RCUT1E=',rcut1e
1756 cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1757 cfmc write(iout,*)'RCUT2S=',rcut2s
1758 cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1759 cfmc write(iout,*)'RCUT2E=',rcut2e
1760 cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1761 cfmc write(iout,*)'DPERT1=',d_pert1
1762 cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1763 cfmc write(iout,*)'DPERT1A=',d_pert1a
1764 cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1765 cfmc write(iout,*)'DPERT2=',d_pert2
1766 cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1767 cfmc write(iout,*)'DPERT2A=',d_pert2a
1768 cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1769 cfmc write(iout,*)'DPERT2B=',d_pert2b
1770 cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1771 cfmc write(iout,*)'DPERT2C=',d_pert2c
1772 cfmc d_pert1=deg2rad*d_pert1
1773 cfmc d_pert1a=deg2rad*d_pert1a
1774 cfmc d_pert2=deg2rad*d_pert2
1775 cfmc d_pert2a=deg2rad*d_pert2a
1776 cfmc d_pert2b=deg2rad*d_pert2b
1777 cfmc d_pert2c=deg2rad*d_pert2c
1778 cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1779 cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
1780 cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1781 cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
1782 cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1783 cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1784 cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1785 cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1786 cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1787 cfmc write(iout,*)'RCUTINI=',rcutini
1788 cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
1789 cfmc write(iout,*)'GRAT=',grat
1790 cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1791 cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
1795 c----------------------------------------------------------------------------
1797 implicit real*8 (a-h,o-z)
1798 include 'DIMENSIONS'
1799 include 'COMMON.MCM'
1800 include 'COMMON.MCE'
1801 include 'COMMON.IOUNITS'
1803 character*320 mcmcard
1804 call card_concat(mcmcard)
1805 call readi(mcmcard,'MAXACC',maxacc,100)
1806 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1807 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1808 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1809 call readi(mcmcard,'MAXREPM',maxrepm,200)
1810 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1811 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1812 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1813 call reada(mcmcard,'E_UP',e_up,5.0D0)
1814 call reada(mcmcard,'DELTE',delte,0.1D0)
1815 call readi(mcmcard,'NSWEEP',nsweep,5)
1816 call readi(mcmcard,'NSTEPH',nsteph,0)
1817 call readi(mcmcard,'NSTEPC',nstepc,0)
1818 call reada(mcmcard,'TMIN',tmin,298.0D0)
1819 call reada(mcmcard,'TMAX',tmax,298.0D0)
1820 call readi(mcmcard,'NWINDOW',nwindow,0)
1821 call readi(mcmcard,'PRINT_MC',print_mc,0)
1822 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1823 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1824 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1825 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1826 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1827 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1828 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1829 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1830 if (nwindow.gt.0) then
1831 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1833 winlen(i)=winend(i)-winstart(i)+1
1836 if (tmax.lt.tmin) tmax=tmin
1837 if (tmax.eq.tmin) then
1841 if (nstepc.gt.0 .and. nsteph.gt.0) then
1842 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1843 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1845 C Probabilities of different move types
1846 sumpro_type(0)=0.0D0
1847 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1848 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1849 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1850 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1851 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1852 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1853 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1855 print *,'i',i,' sumprotype',sumpro_type(i)
1856 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1857 print *,'i',i,' sumprotype',sumpro_type(i)
1861 c----------------------------------------------------------------------------
1862 subroutine read_minim
1863 implicit real*8 (a-h,o-z)
1864 include 'DIMENSIONS'
1865 include 'COMMON.MINIM'
1866 include 'COMMON.IOUNITS'
1868 character*320 minimcard
1869 call card_concat(minimcard)
1870 call readi(minimcard,'MAXMIN',maxmin,2000)
1871 call readi(minimcard,'MAXFUN',maxfun,5000)
1872 call readi(minimcard,'MINMIN',minmin,maxmin)
1873 call readi(minimcard,'MINFUN',minfun,maxmin)
1874 call reada(minimcard,'TOLF',tolf,1.0D-2)
1875 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1876 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1877 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1878 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1879 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1880 & 'Options in energy minimization:'
1881 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1882 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1883 & 'MinMin:',MinMin,' MinFun:',MinFun,
1884 & ' TolF:',TolF,' RTolF:',RTolF
1887 c----------------------------------------------------------------------------
1888 subroutine read_angles(kanal,*)
1889 implicit real*8 (a-h,o-z)
1890 include 'DIMENSIONS'
1891 include 'COMMON.GEO'
1892 include 'COMMON.VAR'
1893 include 'COMMON.CHAIN'
1894 include 'COMMON.IOUNITS'
1895 include 'COMMON.CONTROL'
1896 c Read angles from input
1898 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1899 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1900 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1901 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1904 c 9/7/01 avoid 180 deg valence angle
1905 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1907 theta(i)=deg2rad*theta(i)
1908 phi(i)=deg2rad*phi(i)
1909 alph(i)=deg2rad*alph(i)
1910 omeg(i)=deg2rad*omeg(i)
1915 c----------------------------------------------------------------------------
1916 subroutine reada(rekord,lancuch,wartosc,default)
1918 character*(*) rekord,lancuch
1919 double precision wartosc,default
1922 iread=index(rekord,lancuch)
1923 if (iread.eq.0) then
1927 iread=iread+ilen(lancuch)+1
1928 read (rekord(iread:),*,err=10,end=10) wartosc
1933 c----------------------------------------------------------------------------
1934 subroutine readi(rekord,lancuch,wartosc,default)
1936 character*(*) rekord,lancuch
1937 integer wartosc,default
1940 iread=index(rekord,lancuch)
1941 if (iread.eq.0) then
1945 iread=iread+ilen(lancuch)+1
1946 read (rekord(iread:),*,err=10,end=10) wartosc
1951 c----------------------------------------------------------------------------
1952 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1955 integer tablica(dim),default
1956 character*(*) rekord,lancuch
1963 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1964 if (iread.eq.0) return
1965 iread=iread+ilen(lancuch)+1
1966 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1969 c----------------------------------------------------------------------------
1970 subroutine multreada(rekord,lancuch,tablica,dim,default)
1973 double precision tablica(dim),default
1974 character*(*) rekord,lancuch
1981 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1982 if (iread.eq.0) return
1983 iread=iread+ilen(lancuch)+1
1984 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1987 c----------------------------------------------------------------------------
1988 subroutine openunits
1989 implicit real*8 (a-h,o-z)
1990 include 'DIMENSIONS'
1993 character*16 form,nodename
1996 include 'COMMON.SETUP'
1997 include 'COMMON.IOUNITS'
1999 include 'COMMON.CONTROL'
2000 integer lenpre,lenpot,ilen,lentmp
2002 character*3 out1file_text,ucase
2005 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
2006 call getenv_loc("PREFIX",prefix)
2008 call getenv_loc("POT",pot)
2009 call getenv_loc("DIRTMP",tmpdir)
2010 call getenv_loc("CURDIR",curdir)
2011 call getenv_loc("OUT1FILE",out1file_text)
2012 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
2013 out1file_text=ucase(out1file_text)
2014 if (out1file_text(1:1).eq."Y") then
2017 out1file=fg_rank.gt.0
2022 if (lentmp.gt.0) then
2023 write (*,'(80(1h!))')
2024 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
2025 write (*,'(80(1h!))')
2026 write (*,*)"All output files will be on node /tmp directory."
2028 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
2029 if (me.eq.king) then
2030 write (*,*) "The master node is ",nodename
2031 else if (fg_rank.eq.0) then
2032 write (*,*) "I am the CG slave node ",nodename
2034 write (*,*) "I am the FG slave node ",nodename
2037 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
2038 lenpre = lentmp+lenpre+1
2040 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
2041 C Get the names and open the input files
2042 #if defined(WINIFL) || defined(WINPGI)
2043 open(1,file=pref_orig(:ilen(pref_orig))//
2044 & '.inp',status='old',readonly,shared)
2045 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2046 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2047 C Get parameter filenames and open the parameter files.
2048 call getenv_loc('BONDPAR',bondname)
2049 open (ibond,file=bondname,status='old',readonly,shared)
2050 call getenv_loc('THETPAR',thetname)
2051 open (ithep,file=thetname,status='old',readonly,shared)
2053 call getenv_loc('THETPARPDB',thetname_pdb)
2054 open (ithep_pdb,file=thetname_pdb,status='old',readonly,shared)
2056 call getenv_loc('ROTPAR',rotname)
2057 open (irotam,file=rotname,status='old',readonly,shared)
2059 call getenv_loc('ROTPARPDB',rotname_pdb)
2060 open (irotam_pdb,file=rotname_pdb,status='old',readonly,shared)
2062 call getenv_loc('TORPAR',torname)
2063 open (itorp,file=torname,status='old',readonly,shared)
2064 call getenv_loc('TORDPAR',tordname)
2065 open (itordp,file=tordname,status='old',readonly,shared)
2066 call getenv_loc('FOURIER',fouriername)
2067 open (ifourier,file=fouriername,status='old',readonly,shared)
2068 call getenv_loc('ELEPAR',elename)
2069 open (ielep,file=elename,status='old',readonly,shared)
2070 call getenv_loc('SIDEPAR',sidename)
2071 open (isidep,file=sidename,status='old',readonly,shared)
2072 #elif (defined CRAY) || (defined AIX)
2073 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2075 c print *,"Processor",myrank," opened file 1"
2076 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2077 c print *,"Processor",myrank," opened file 9"
2078 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2079 C Get parameter filenames and open the parameter files.
2080 call getenv_loc('BONDPAR',bondname)
2081 open (ibond,file=bondname,status='old',action='read')
2082 c print *,"Processor",myrank," opened file IBOND"
2083 call getenv_loc('THETPAR',thetname)
2084 open (ithep,file=thetname,status='old',action='read')
2085 c print *,"Processor",myrank," opened file ITHEP"
2087 call getenv_loc('THETPARPDB',thetname_pdb)
2088 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2090 call getenv_loc('ROTPAR',rotname)
2091 open (irotam,file=rotname,status='old',action='read')
2092 c print *,"Processor",myrank," opened file IROTAM"
2094 call getenv_loc('ROTPARPDB',rotname_pdb)
2095 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2097 call getenv_loc('TORPAR',torname)
2098 open (itorp,file=torname,status='old',action='read')
2099 c print *,"Processor",myrank," opened file ITORP"
2100 call getenv_loc('TORDPAR',tordname)
2101 open (itordp,file=tordname,status='old',action='read')
2102 c print *,"Processor",myrank," opened file ITORDP"
2103 call getenv_loc('SCCORPAR',sccorname)
2104 open (isccor,file=sccorname,status='old',action='read')
2105 c print *,"Processor",myrank," opened file ISCCOR"
2106 call getenv_loc('FOURIER',fouriername)
2107 open (ifourier,file=fouriername,status='old',action='read')
2108 c print *,"Processor",myrank," opened file IFOURIER"
2109 call getenv_loc('ELEPAR',elename)
2110 open (ielep,file=elename,status='old',action='read')
2111 c print *,"Processor",myrank," opened file IELEP"
2112 call getenv_loc('SIDEPAR',sidename)
2113 open (isidep,file=sidename,status='old',action='read')
2114 c print *,"Processor",myrank," opened file ISIDEP"
2115 c print *,"Processor",myrank," opened parameter files"
2117 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2118 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2119 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2120 C Get parameter filenames and open the parameter files.
2121 call getenv_loc('BONDPAR',bondname)
2122 open (ibond,file=bondname,status='old')
2123 call getenv_loc('THETPAR',thetname)
2124 open (ithep,file=thetname,status='old')
2126 call getenv_loc('THETPARPDB',thetname_pdb)
2127 open (ithep_pdb,file=thetname_pdb,status='old')
2129 call getenv_loc('ROTPAR',rotname)
2130 open (irotam,file=rotname,status='old')
2132 call getenv_loc('ROTPARPDB',rotname_pdb)
2133 open (irotam_pdb,file=rotname_pdb,status='old')
2135 call getenv_loc('TORPAR',torname)
2136 open (itorp,file=torname,status='old')
2137 call getenv_loc('TORDPAR',tordname)
2138 open (itordp,file=tordname,status='old')
2139 call getenv_loc('SCCORPAR',sccorname)
2140 open (isccor,file=sccorname,status='old')
2141 call getenv_loc('FOURIER',fouriername)
2142 open (ifourier,file=fouriername,status='old')
2143 call getenv_loc('ELEPAR',elename)
2144 open (ielep,file=elename,status='old')
2145 call getenv_loc('SIDEPAR',sidename)
2146 open (isidep,file=sidename,status='old')
2148 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2150 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2151 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2152 C Get parameter filenames and open the parameter files.
2153 call getenv_loc('BONDPAR',bondname)
2154 open (ibond,file=bondname,status='old',action='read')
2155 call getenv_loc('THETPAR',thetname)
2156 open (ithep,file=thetname,status='old',action='read')
2158 call getenv_loc('THETPARPDB',thetname_pdb)
2159 print *,"thetname_pdb ",thetname_pdb
2160 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2161 print *,ithep_pdb," opened"
2163 call getenv_loc('ROTPAR',rotname)
2164 open (irotam,file=rotname,status='old',action='read')
2166 call getenv_loc('ROTPARPDB',rotname_pdb)
2167 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2169 call getenv_loc('TORPAR',torname)
2170 open (itorp,file=torname,status='old',action='read')
2171 call getenv_loc('TORDPAR',tordname)
2172 open (itordp,file=tordname,status='old',action='read')
2173 call getenv_loc('SCCORPAR',sccorname)
2174 open (isccor,file=sccorname,status='old',action='read')
2175 call getenv_loc('FOURIER',fouriername)
2176 open (ifourier,file=fouriername,status='old',action='read')
2177 call getenv_loc('ELEPAR',elename)
2178 open (ielep,file=elename,status='old',action='read')
2179 call getenv_loc('SIDEPAR',sidename)
2180 open (isidep,file=sidename,status='old',action='read')
2184 C 8/9/01 In the newest version SCp interaction constants are read from a file
2185 C Use -DOLDSCP to use hard-coded constants instead.
2187 call getenv_loc('SCPPAR',scpname)
2188 #if defined(WINIFL) || defined(WINPGI)
2189 open (iscpp,file=scpname,status='old',readonly,shared)
2190 #elif (defined CRAY) || (defined AIX)
2191 open (iscpp,file=scpname,status='old',action='read')
2193 open (iscpp,file=scpname,status='old')
2195 open (iscpp,file=scpname,status='old',action='read')
2198 call getenv_loc('PATTERN',patname)
2199 #if defined(WINIFL) || defined(WINPGI)
2200 open (icbase,file=patname,status='old',readonly,shared)
2201 #elif (defined CRAY) || (defined AIX)
2202 open (icbase,file=patname,status='old',action='read')
2204 open (icbase,file=patname,status='old')
2206 open (icbase,file=patname,status='old',action='read')
2209 C Open output file only for CG processes
2210 c print *,"Processor",myrank," fg_rank",fg_rank
2211 if (fg_rank.eq.0) then
2213 if (nodes.eq.1) then
2216 npos = dlog10(dfloat(nodes-1))+1
2218 if (npos.lt.3) npos=3
2219 write (liczba,'(i1)') npos
2220 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2222 write (liczba,form) me
2223 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2224 & liczba(:ilen(liczba))
2225 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2227 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2229 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2230 & liczba(:ilen(liczba))//'.mol2'
2231 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2232 & liczba(:ilen(liczba))//'.stat'
2234 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2235 & //liczba(:ilen(liczba))//'.stat')
2236 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2239 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2240 & liczba(:ilen(liczba))//'.const'
2245 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2246 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2247 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2248 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2249 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2251 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2253 rest2name=prefix(:ilen(prefix))//'.rst'
2255 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2258 #if defined(AIX) || defined(PGI)
2259 if (me.eq.king .or. .not. out1file)
2260 & open(iout,file=outname,status='unknown')
2263 if (fg_rank.gt.0) then
2264 write (liczba,'(i3.3)') myrank/nfgtasks
2265 write (ll,'(bz,i3.3)') fg_rank
2266 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2272 open(igeom,file=intname,status='unknown',position='append')
2273 open(ipdb,file=pdbname,status='unknown')
2274 open(imol2,file=mol2name,status='unknown')
2275 open(istat,file=statname,status='unknown',position='append')
2277 c1out open(iout,file=outname,status='unknown')
2280 if (me.eq.king .or. .not.out1file)
2281 & open(iout,file=outname,status='unknown')
2284 if (fg_rank.gt.0) then
2285 print "Processor",fg_rank," opening output file"
2286 write (liczba,'(i3.3)') myrank/nfgtasks
2287 write (ll,'(bz,i3.3)') fg_rank
2288 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2294 open(igeom,file=intname,status='unknown',access='append')
2295 open(ipdb,file=pdbname,status='unknown')
2296 open(imol2,file=mol2name,status='unknown')
2297 open(istat,file=statname,status='unknown',access='append')
2299 c1out open(iout,file=outname,status='unknown')
2302 csa csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2303 csa csa_seed=prefix(:lenpre)//'.CSA.seed'
2304 csa csa_history=prefix(:lenpre)//'.CSA.history'
2305 csa csa_bank=prefix(:lenpre)//'.CSA.bank'
2306 csa csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2307 csa csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2308 csa csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2309 csac!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2310 csa csa_int=prefix(:lenpre)//'.int'
2311 csa csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2312 csa csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2313 csa csa_in=prefix(:lenpre)//'.CSA.in'
2314 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2317 write (iout,'(80(1h-))')
2318 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2319 write (iout,'(80(1h-))')
2320 write (iout,*) "Input file : ",
2321 & pref_orig(:ilen(pref_orig))//'.inp'
2322 write (iout,*) "Output file : ",
2323 & outname(:ilen(outname))
2325 write (iout,*) "Sidechain potential file : ",
2326 & sidename(:ilen(sidename))
2328 write (iout,*) "SCp potential file : ",
2329 & scpname(:ilen(scpname))
2331 write (iout,*) "Electrostatic potential file : ",
2332 & elename(:ilen(elename))
2333 write (iout,*) "Cumulant coefficient file : ",
2334 & fouriername(:ilen(fouriername))
2335 write (iout,*) "Torsional parameter file : ",
2336 & torname(:ilen(torname))
2337 write (iout,*) "Double torsional parameter file : ",
2338 & tordname(:ilen(tordname))
2339 write (iout,*) "SCCOR parameter file : ",
2340 & sccorname(:ilen(sccorname))
2341 write (iout,*) "Bond & inertia constant file : ",
2342 & bondname(:ilen(bondname))
2343 write (iout,*) "Bending parameter file : ",
2344 & thetname(:ilen(thetname))
2345 write (iout,*) "Rotamer parameter file : ",
2346 & rotname(:ilen(rotname))
2347 write (iout,*) "Threading database : ",
2348 & patname(:ilen(patname))
2350 &write (iout,*)" DIRTMP : ",
2352 write (iout,'(80(1h-))')
2356 c----------------------------------------------------------------------------
2357 subroutine card_concat(card)
2358 implicit real*8 (a-h,o-z)
2359 include 'DIMENSIONS'
2360 include 'COMMON.IOUNITS'
2362 character*80 karta,ucase
2364 read (inp,'(a)') karta
2367 do while (karta(80:80).eq.'&')
2368 card=card(:ilen(card)+1)//karta(:79)
2369 read (inp,'(a)') karta
2372 card=card(:ilen(card)+1)//karta
2375 c----------------------------------------------------------------------------------
2377 implicit real*8 (a-h,o-z)
2378 include 'DIMENSIONS'
2379 include 'COMMON.CHAIN'
2380 include 'COMMON.IOUNITS'
2382 open(irest2,file=rest2name,status='unknown')
2383 read(irest2,*) totT,EK,potE,totE,t_bath
2385 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2388 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2391 read (irest2,*) iset
2396 c---------------------------------------------------------------------------------
2397 subroutine read_fragments
2398 implicit real*8 (a-h,o-z)
2399 include 'DIMENSIONS'
2403 include 'COMMON.SETUP'
2404 include 'COMMON.CHAIN'
2405 include 'COMMON.IOUNITS'
2407 include 'COMMON.CONTROL'
2408 read(inp,*) nset,nfrag,npair,nfrag_back
2409 if(me.eq.king.or..not.out1file)
2410 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2411 & " nfrag_back",nfrag_back
2413 read(inp,*) mset(iset)
2415 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2417 if(me.eq.king.or..not.out1file)
2418 & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2419 & ifrag(2,i,iset), qinfrag(i,iset)
2422 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2424 if(me.eq.king.or..not.out1file)
2425 & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2426 & ipair(2,i,iset), qinpair(i,iset)
2429 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2430 & wfrag_back(3,i,iset),
2431 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2432 if(me.eq.king.or..not.out1file)
2433 & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2434 & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2439 c-------------------------------------------------------------------------------
2440 subroutine read_dist_constr
2441 implicit real*8 (a-h,o-z)
2442 include 'DIMENSIONS'
2446 include 'COMMON.SETUP'
2447 include 'COMMON.CONTROL'
2448 include 'COMMON.CHAIN'
2449 include 'COMMON.IOUNITS'
2450 include 'COMMON.SBRIDGE'
2451 integer ifrag_(2,100),ipair_(2,100)
2452 double precision wfrag_(100),wpair_(100)
2453 character*500 controlcard
2454 c write (iout,*) "Calling read_dist_constr"
2455 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2457 call card_concat(controlcard)
2458 call readi(controlcard,"NFRAG",nfrag_,0)
2459 call readi(controlcard,"NPAIR",npair_,0)
2460 call readi(controlcard,"NDIST",ndist_,0)
2461 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2462 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2463 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2464 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2465 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2466 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2467 c write (iout,*) "IFRAG"
2469 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2471 c write (iout,*) "IPAIR"
2473 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2475 if (.not.refstr .and. nfrag.gt.0) then
2477 & "ERROR: no reference structure to compute distance restraints"
2479 & "Restraints must be specified explicitly (NDIST=number)"
2482 if (nfrag.lt.2 .and. npair.gt.0) then
2483 write (iout,*) "ERROR: Less than 2 fragments specified",
2484 & " but distance restraints between pairs requested"
2489 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2490 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2491 & ifrag_(2,i)=nstart_sup+nsup-1
2492 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2494 if (wfrag_(i).gt.0.0d0) then
2495 do j=ifrag_(1,i),ifrag_(2,i)-1
2496 do k=j+1,ifrag_(2,i)
2497 c write (iout,*) "j",j," k",k
2499 if (constr_dist.eq.1) then
2504 forcon(nhpb)=wfrag_(i)
2505 else if (constr_dist.eq.2) then
2506 if (ddjk.le.dist_cut) then
2511 forcon(nhpb)=wfrag_(i)
2518 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2521 if (.not.out1file .or. me.eq.king)
2522 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2523 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2525 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2526 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2533 if (wpair_(i).gt.0.0d0) then
2541 do j=ifrag_(1,ii),ifrag_(2,ii)
2542 do k=ifrag_(1,jj),ifrag_(2,jj)
2546 forcon(nhpb)=wpair_(i)
2547 dhpb(nhpb)=dist(j,k)
2549 if (.not.out1file .or. me.eq.king)
2550 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2551 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2553 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2554 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2561 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2562 & ibecarb(i),forcon(nhpb+1)
2563 if (forcon(nhpb+1).gt.0.0d0) then
2565 if (ibecarb(i).gt.0) then
2566 ihpb(i)=ihpb(i)+nres
2567 jhpb(i)=jhpb(i)+nres
2569 if (dhpb(nhpb).eq.0.0d0)
2570 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2574 if (.not.out1file .or. me.eq.king) then
2577 write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ",
2578 & i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i)
2586 c-------------------------------------------------------------------------------
2588 subroutine read_constr_homology
2590 include 'DIMENSIONS'
2594 include 'COMMON.SETUP'
2595 include 'COMMON.CONTROL'
2596 include 'COMMON.CHAIN'
2597 include 'COMMON.IOUNITS'
2601 character*24 model_ki_dist, model_ki_angle
2602 character*500 controlcard
2603 integer ki, i, j, k, l
2606 call card_concat(controlcard)
2607 call reada(controlcard,"HOMOL_DIST",waga_dist,1.0)
2608 call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0)
2610 do ki=1,constr_homology
2611 if (constr_homology.ge.1) then
2613 write(kic2,'(i2)') ki
2614 c write(iout,*) "TEST KICA, HOMOL", kic2
2615 if (ki.le.9) kic2="0"//kic2(2:2)
2616 c write(iout,*) "TEST KICA2, HOMOL", kic2
2618 model_ki_dist="model"//kic2//".dist"
2619 model_ki_angle="model"//kic2//".angle"
2620 c write(iout,*) model_ki_dist, model_ki_angle
2621 open (1400+ki,file=model_ki_dist,status='old')
2622 open (1401+ki,file=model_ki_angle,status='old')
2624 do irec=1,99999 !petla do czytania wiezow na odleglosc
2625 read (1400+ki,*,end=1401) i, j, odl(i,j,ki),sigma_odl(i,j,ki)
2629 do irec=1,99999 !petla do czytania wiezow na katach torsyjnych
2630 read (1401+ki,*,end=1402) i, j, k,l,dih(i,ki),sigma_dih(i,ki)
2632 c dih(i,ki)=dih(i,ki)
2638 c write(iout,*) "TEST CZYTANIA1",odl(1,2,1),odl(1,3,1),odl(1,4,1)
2639 c write(iout,*) "TEST CZYTANIA2",dih(1,1),dih(2,1),dih(3,1)
2644 c----------------------------------------------------------------------
2647 subroutine flush(iu)
2652 subroutine flush(iu)
2657 c------------------------------------------------------------------------------
2658 subroutine copy_to_tmp(source)
2659 include "DIMENSIONS"
2660 include "COMMON.IOUNITS"
2661 character*(*) source
2662 character* 256 tmpfile
2666 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
2667 inquire(file=tmpfile,exist=ex)
2669 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
2670 & " to temporary directory..."
2671 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
2672 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
2676 c------------------------------------------------------------------------------
2677 subroutine move_from_tmp(source)
2678 include "DIMENSIONS"
2679 include "COMMON.IOUNITS"
2680 character*(*) source
2683 write (*,*) "Moving ",source(:ilen(source)),
2684 & " from temporary directory to working directory"
2685 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
2686 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
2689 c------------------------------------------------------------------------------
2690 subroutine random_init(seed)
2692 C Initialize random number generator
2694 implicit real*8 (a-h,o-z)
2695 include 'DIMENSIONS'
2701 logical OKRandom, prng_restart
2703 integer iseed_array(4)
2705 include 'COMMON.IOUNITS'
2706 include 'COMMON.TIME1'
2707 include 'COMMON.THREAD'
2708 include 'COMMON.SBRIDGE'
2709 include 'COMMON.CONTROL'
2710 include 'COMMON.MCM'
2711 include 'COMMON.MAP'
2712 include 'COMMON.HEADER'
2713 csa include 'COMMON.CSA'
2714 include 'COMMON.CHAIN'
2715 include 'COMMON.MUCA'
2717 include 'COMMON.FFIELD'
2718 include 'COMMON.SETUP'
2719 iseed=-dint(dabs(seed))
2720 if (iseed.eq.0) then
2721 write (iout,'(/80(1h*)/20x,a/80(1h*))')
2722 & 'Random seed undefined. The program will stop.'
2723 write (*,'(/80(1h*)/20x,a/80(1h*))')
2724 & 'Random seed undefined. The program will stop.'
2726 call mpi_finalize(mpi_comm_world,ierr)
2728 stop 'Bad random seed.'
2731 if (fg_rank.eq.0) then
2735 if(me.eq.king .or. .not. out1file)
2736 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8
2737 write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8
2738 OKRandom = prng_restart(me,iseedi8)
2741 tmp=65536.0d0**(4-i)
2742 iseed_array(i) = dint(seed/tmp)
2743 seed=seed-iseed_array(i)*tmp
2745 if(me.eq.king .or. .not. out1file)
2746 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
2747 & (iseed_array(i),i=1,4)
2748 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
2749 & (iseed_array(i),i=1,4)
2750 OKRandom = prng_restart(me,iseed_array)
2753 r1=ran_number(0.0D0,1.0D0)
2754 if(me.eq.king .or. .not. out1file)
2755 & write (iout,*) 'ran_num',r1
2756 if (r1.lt.0.0d0) OKRandom=.false.
2758 if (.not.OKRandom) then
2759 write (iout,*) 'PRNG IS NOT WORKING!!!'
2760 print *,'PRNG IS NOT WORKING!!!'
2763 call mpi_abort(mpi_comm_world,error_msg,ierr)
2766 write (iout,*) 'too many processors for parallel prng'
2767 write (*,*) 'too many processors for parallel prng'
2775 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)