2 implicit real*8 (a-h,o-z)
8 include 'COMMON.CONTROL'
9 include 'COMMON.SBRIDGE'
10 include 'COMMON.IOUNITS'
11 include 'COMMON.CHAIN'
13 C Read force-field parameters except weights
15 C Read job setup parameters
17 C Read control parameters for energy minimzation if required
18 if (minim) call read_minim
19 C Read MCM control parameters if required
20 if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
21 C Read MD control parameters if reqjuired
22 if (modecalc.eq.12) call read_MDpar
23 C Read MREMD control parameters if required
24 if (modecalc.eq.14) then
28 C Read MUCA control parameters if required
29 if (lmuca) call read_muca
30 C Read CSA control parameters if required (from fort.40 if exists
31 C otherwise from general input file)
32 csa if (modecalc.eq.8) then
33 csa inquire (file="fort.40",exist=file_exist)
34 csa if (.not.file_exist) call csaread
36 cfmc if (modecalc.eq.10) call mcmfread
37 C Read molecule information, molecule geometry, energy-term weights, and
38 C restraints if requested
40 C Print restraint information
42 if (.not. out1file .or. me.eq.king) then
45 &write (iout,'(a,i5,a)') "The following",nhpb-nss,
46 & " distance constraints have been imposed"
48 write (iout,'(3i6,i2,3f10.5)') i-nss,ihpb(i),jhpb(i),
49 & ibecarb(i),dhpb(i),dhpb1(i),forcon(i)
54 c print *,"Processor",myrank," leaves READRTNS"
57 C-------------------------------------------------------------------------------
58 subroutine read_control
62 implicit real*8 (a-h,o-z)
66 logical OKRandom, prng_restart
69 include 'COMMON.IOUNITS'
70 include 'COMMON.TIME1'
71 include 'COMMON.THREAD'
72 include 'COMMON.SBRIDGE'
73 include 'COMMON.CONTROL'
76 include 'COMMON.HEADER'
77 csa include 'COMMON.CSA'
78 include 'COMMON.CHAIN'
81 include 'COMMON.FFIELD'
82 include 'COMMON.SETUP'
83 COMMON /MACHSW/ KDIAG,ICORFL,IXDR
84 character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
86 character*320 controlcard
91 read (INP,'(a)') titel
92 call card_concat(controlcard)
93 c out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
94 c print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
95 call reada(controlcard,'SEED',seed,0.0D0)
96 call random_init(seed)
97 C Set up the time limit (caution! The time must be input in minutes!)
98 read_cart=index(controlcard,'READ_CART').gt.0
99 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
100 call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
101 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
102 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
103 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
104 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
105 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
106 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
107 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
108 call reada(controlcard,'DRMS',drms,0.1D0)
109 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
110 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
111 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
112 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
113 write (iout,'(a,f10.1)')'DRMS = ',drms
114 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
115 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
117 call readi(controlcard,'NZ_START',nz_start,0)
118 call readi(controlcard,'NZ_END',nz_end,0)
119 call readi(controlcard,'IZ_SC',iz_sc,0)
121 safety = 60.0d0*safety
124 call reada(controlcard,"T_BATH",t_bath,300.0d0)
125 minim=(index(controlcard,'MINIMIZE').gt.0)
126 dccart=(index(controlcard,'CART').gt.0)
127 overlapsc=(index(controlcard,'OVERLAP').gt.0)
128 overlapsc=.not.overlapsc
129 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
130 searchsc=.not.searchsc
131 sideadd=(index(controlcard,'SIDEADD').gt.0)
132 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
133 outpdb=(index(controlcard,'PDBOUT').gt.0)
134 outx=(index(controlcard,'XOUT').gt.0)
135 outmol2=(index(controlcard,'MOL2OUT').gt.0)
136 pdbref=(index(controlcard,'PDBREF').gt.0)
137 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
138 indpdb=index(controlcard,'PDBSTART')
139 extconf=(index(controlcard,'EXTCONF').gt.0)
140 call readi(controlcard,'IPRINT',iprint,0)
141 call readi(controlcard,'MAXGEN',maxgen,10000)
142 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
143 call readi(controlcard,"KDIAG",kdiag,0)
144 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
145 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
146 & write (iout,*) "RESCALE_MODE",rescale_mode
147 split_ene=index(controlcard,'SPLIT_ENE').gt.0
148 if (index(controlcard,'REGULAR').gt.0.0D0) then
149 call reada(controlcard,'WEIDIS',weidis,0.1D0)
153 if (index(controlcard,'CHECKGRAD').gt.0) then
155 if (index(controlcard,'CART').gt.0) then
157 elseif (index(controlcard,'CARINT').gt.0) then
162 elseif (index(controlcard,'THREAD').gt.0) then
164 call readi(controlcard,'THREAD',nthread,0)
165 if (nthread.gt.0) then
166 call reada(controlcard,'WEIDIS',weidis,0.1D0)
169 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
170 stop 'Error termination in Read_Control.'
172 else if (index(controlcard,'MCMA').gt.0) then
174 else if (index(controlcard,'MCEE').gt.0) then
176 else if (index(controlcard,'MULTCONF').gt.0) then
178 else if (index(controlcard,'MAP').gt.0) then
180 call readi(controlcard,'MAP',nmap,0)
181 else if (index(controlcard,'CSA').gt.0) then
182 write(*,*) "CSA not supported in this version"
185 crc else if (index(controlcard,'ZSCORE').gt.0) then
187 crc ZSCORE is rm from UNRES, modecalc=9 is available
190 cfcm else if (index(controlcard,'MCMF').gt.0) then
192 else if (index(controlcard,'SOFTREG').gt.0) then
194 else if (index(controlcard,'CHECK_BOND').gt.0) then
196 else if (index(controlcard,'TEST').gt.0) then
198 else if (index(controlcard,'MD').gt.0) then
200 else if (index(controlcard,'RE ').gt.0) then
204 lmuca=index(controlcard,'MUCA').gt.0
205 call readi(controlcard,'MUCADYN',mucadyn,0)
206 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
207 if (lmuca .and. (me.eq.king .or. .not.out1file ))
209 write (iout,*) 'MUCADYN=',mucadyn
210 write (iout,*) 'MUCASMOOTH=',muca_smooth
213 iscode=index(controlcard,'ONE_LETTER')
214 indphi=index(controlcard,'PHI')
215 indback=index(controlcard,'BACK')
216 iranconf=index(controlcard,'RAND_CONF')
217 i2ndstr=index(controlcard,'USE_SEC_PRED')
218 gradout=index(controlcard,'GRADOUT').gt.0
219 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
221 if(me.eq.king.or..not.out1file)
222 & write (iout,'(2a)') diagmeth(kdiag),
223 & ' routine used to diagonalize matrices.'
226 c--------------------------------------------------------------------------
227 subroutine read_REMDpar
231 implicit real*8 (a-h,o-z)
233 include 'COMMON.IOUNITS'
234 include 'COMMON.TIME1'
237 include 'COMMON.LANGEVIN'
239 include 'COMMON.LANGEVIN.lang0'
241 include 'COMMON.INTERACT'
242 include 'COMMON.NAMES'
244 include 'COMMON.REMD'
245 include 'COMMON.CONTROL'
246 include 'COMMON.SETUP'
248 character*320 controlcard
249 character*3200 controlcard1
250 integer iremd_m_total
252 if(me.eq.king.or..not.out1file)
253 & write (iout,*) "REMD setup"
255 call card_concat(controlcard)
256 call readi(controlcard,"NREP",nrep,3)
257 call readi(controlcard,"NSTEX",nstex,1000)
258 call reada(controlcard,"RETMIN",retmin,10.0d0)
259 call reada(controlcard,"RETMAX",retmax,1000.0d0)
260 mremdsync=(index(controlcard,'SYNC').gt.0)
261 call readi(controlcard,"NSYN",i_sync_step,100)
262 restart1file=(index(controlcard,'REST1FILE').gt.0)
263 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
264 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
265 if(max_cache_traj_use.gt.max_cache_traj)
266 & max_cache_traj_use=max_cache_traj
267 if(me.eq.king.or..not.out1file) then
268 cd if (traj1file) then
269 crc caching is in testing - NTWX is not ignored
270 cd write (iout,*) "NTWX value is ignored"
271 cd write (iout,*) " trajectory is stored to one file by master"
272 cd write (iout,*) " before exchange at NSTEX intervals"
274 write (iout,*) "NREP= ",nrep
275 write (iout,*) "NSTEX= ",nstex
276 write (iout,*) "SYNC= ",mremdsync
277 write (iout,*) "NSYN= ",i_sync_step
278 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
281 t_exchange_only=(index(controlcard,'TONLY').gt.0)
282 call readi(controlcard,"HREMD",hremd,0)
283 if((me.eq.king.or..not.out1file).and.hremd.gt.0) then
284 write (iout,*) "Hamiltonian REMD with ",hremd," sets of weights"
286 if(usampl.and.hremd.gt.0) then
288 & "========== ERROR: USAMPL and HREMD cannot be used together"
290 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
297 if (index(controlcard,'TLIST').gt.0) then
299 call card_concat(controlcard1)
300 read(controlcard1,*) (remd_t(i),i=1,nrep)
301 if(me.eq.king.or..not.out1file)
302 & write (iout,*)'tlist',(remd_t(i),i=1,nrep)
305 if (index(controlcard,'MLIST').gt.0) then
307 call card_concat(controlcard1)
308 read(controlcard1,*) (remd_m(i),i=1,nrep)
309 if(me.eq.king.or..not.out1file) then
310 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
313 iremd_m_total=iremd_m_total+remd_m(i)
316 write (iout,*) 'Total number of replicas ',
317 & iremd_m_total*hremd
319 write (iout,*) 'Total number of replicas ',iremd_m_total
323 if(me.eq.king.or..not.out1file)
324 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
327 c--------------------------------------------------------------------------
328 subroutine read_MDpar
332 implicit real*8 (a-h,o-z)
334 include 'COMMON.IOUNITS'
335 include 'COMMON.TIME1'
338 include 'COMMON.LANGEVIN'
340 include 'COMMON.LANGEVIN.lang0'
342 include 'COMMON.INTERACT'
343 include 'COMMON.NAMES'
345 include 'COMMON.SETUP'
346 include 'COMMON.CONTROL'
347 include 'COMMON.SPLITELE'
349 character*320 controlcard
351 call card_concat(controlcard)
352 call readi(controlcard,"NSTEP",n_timestep,1000000)
353 call readi(controlcard,"NTWE",ntwe,100)
354 call readi(controlcard,"NTWX",ntwx,1000)
355 call reada(controlcard,"DT",d_time,1.0d-1)
356 call reada(controlcard,"DVMAX",dvmax,2.0d1)
357 call reada(controlcard,"DAMAX",damax,1.0d1)
358 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
359 call readi(controlcard,"LANG",lang,0)
360 RESPA = index(controlcard,"RESPA") .gt. 0
361 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
362 ntime_split0=ntime_split
363 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
364 ntime_split0=ntime_split
365 call reada(controlcard,"R_CUT",r_cut,2.0d0)
366 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
367 rest = index(controlcard,"REST").gt.0
368 tbf = index(controlcard,"TBF").gt.0
369 call readi(controlcard,"HMC",hmc,0)
370 tnp = index(controlcard,"NOSEPOINCARE99").gt.0
371 tnp1 = index(controlcard,"NOSEPOINCARE01").gt.0
372 tnh = index(controlcard,"NOSEHOOVER96").gt.0
373 if (RESPA.and.tnh)then
374 xiresp = index(controlcard,"XIRESP").gt.0
376 call reada(controlcard,"Q_NP",Q_np,0.1d0)
377 usampl = index(controlcard,"USAMPL").gt.0
378 mdpdb = index(controlcard,"MDPDB").gt.0
379 call reada(controlcard,"T_BATH",t_bath,300.0d0)
380 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
381 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
382 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
383 if (count_reset_moment.eq.0) count_reset_moment=1000000000
384 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
385 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
386 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
387 if (count_reset_vel.eq.0) count_reset_vel=1000000000
388 large = index(controlcard,"LARGE").gt.0
389 print_compon = index(controlcard,"PRINT_COMPON").gt.0
390 rattle = index(controlcard,"RATTLE").gt.0
391 preminim = index(controlcard,"PREMINIM").gt.0
393 dccart=(index(controlcard,'CART').gt.0)
396 c if performing umbrella sampling, fragments constrained are read from the fragment file
402 if(me.eq.king.or..not.out1file) then
404 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
406 write (iout,'(a)') "The units are:"
407 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
408 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
409 & " acceleration: angstrom/(48.9 fs)**2"
410 write (iout,'(a)') "energy: kcal/mol, temperature: K"
412 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
413 write (iout,'(a60,f10.5,a)')
414 & "Initial time step of numerical integration:",d_time,
416 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
418 write (iout,'(2a,i4,a)')
419 & "A-MTS algorithm used; initial time step for fast-varying",
420 & " short-range forces split into",ntime_split," steps."
421 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
422 & r_cut," lambda",rlamb
424 write (iout,'(2a,f10.5)')
425 & "Maximum acceleration threshold to reduce the time step",
426 & "/increase split number:",damax
427 write (iout,'(2a,f10.5)')
428 & "Maximum predicted energy drift to reduce the timestep",
429 & "/increase split number:",edriftmax
430 write (iout,'(a60,f10.5)')
431 & "Maximum velocity threshold to reduce velocities:",dvmax
432 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
433 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
434 if (rattle) write (iout,'(a60)')
435 & "Rattle algorithm used to constrain the virtual bonds"
436 if (preminim .or. iranconf.gt.0) then
438 & "Initial structure will be energy-minimized"
443 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
444 call reada(controlcard,"RWAT",rwat,1.4d0)
445 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
446 surfarea=index(controlcard,"SURFAREA").gt.0
447 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
448 if(me.eq.king.or..not.out1file)then
449 write (iout,'(/a,$)') "Langevin dynamics calculation"
452 & " with direct integration of Langevin equations"
453 else if (lang.eq.2) then
454 write (iout,'(a/)') " with TINKER stochasic MD integrator"
455 else if (lang.eq.3) then
456 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
457 else if (lang.eq.4) then
458 write (iout,'(a/)') " in overdamped mode"
460 write (iout,'(//a,i5)')
461 & "=========== ERROR: Unknown Langevin dynamics mode:",lang
464 write (iout,'(a60,f10.5)') "Temperature:",t_bath
465 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
466 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
467 write (iout,'(a60,f10.5)')
468 & "Scaling factor of the friction forces:",scal_fric
469 if (surfarea) write (iout,'(2a,i10,a)')
470 & "Friction coefficients will be scaled by solvent-accessible",
471 & " surface area every",reset_fricmat," steps."
473 c Calculate friction coefficients and bounds of stochastic forces
474 eta=6*pi*cPoise*etawat
475 if(me.eq.king.or..not.out1file)
476 & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
478 gamp=scal_fric*(pstok+rwat)*eta
479 stdfp=dsqrt(2*Rb*t_bath/d_time)
481 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
482 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
484 if(me.eq.king.or..not.out1file)then
485 write (iout,'(/2a/)')
486 & "Radii of site types and friction coefficients and std's of",
487 & " stochastic forces of fully exposed sites"
488 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
490 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
491 & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
495 if(me.eq.king.or..not.out1file)then
496 write (iout,'(a)') "Berendsen bath calculation"
497 write (iout,'(a60,f10.5)') "Temperature:",t_bath
498 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
500 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
501 & count_reset_moment," steps"
503 & write (iout,'(a,i10,a)')
504 & "Velocities will be reset at random every",count_reset_vel,
507 else if (tnp .or. tnp1 .or. tnh) then
508 if (tnp .or. tnp1) then
509 write (iout,'(a)') "Nose-Poincare bath calculation"
510 if (tnp) write (iout,'(a)')
511 & "J.Comput.Phys. 151 114 (1999) S.D.Bond B.J.Leimkuhler B.B.Laird"
512 if (tnp1) write (iout,'(a)') "JPSJ 70 75 (2001) S. Nose"
514 write (iout,'(a)') "Nose-Hoover bath calculation"
515 write (iout,'(a)') "Mol.Phys. 87 1117 (1996) Martyna et al."
525 WDTI(i) = 1.0*d_time/nresn
532 write (iout,'(a)') "NVT-XI-RESPA algorithm"
534 write (iout,'(a)') "NVT-XO-RESPA algorithm"
537 WDTIi(i) = 1.0*d_time/nresn/ntime_split
545 write (iout,'(a60,f10.5)') "Temperature:",t_bath
546 write (iout,'(a60,f10.5)') "Q =",Q_np
548 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
549 & count_reset_moment," steps"
551 & write (iout,'(a,i10,a)')
552 & "Velocities will be reset at random every",count_reset_vel,
555 else if (hmc.gt.0) then
556 write (iout,'(a)') "Hybrid Monte Carlo calculation"
557 write (iout,'(a60,f10.5)') "Temperature:",t_bath
558 write (iout,'(a60,i10)')
559 & "Number of MD steps between Metropolis tests:",hmc
562 if(me.eq.king.or..not.out1file)
563 & write (iout,'(a31)') "Microcanonical mode calculation"
565 if(me.eq.king.or..not.out1file)then
566 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
568 write(iout,*) "MD running with constraints."
569 write(iout,*) "Equilibration time ", eq_time, " mtus."
570 write(iout,*) "Constraining ", nfrag," fragments."
571 write(iout,*) "Length of each fragment, weight and q0:"
573 write (iout,*) "Set of restraints #",iset
575 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
576 & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
578 write(iout,*) "constraints between ", npair, "fragments."
579 write(iout,*) "constraint pairs, weights and q0:"
581 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
582 & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
584 write(iout,*) "angle constraints within ", nfrag_back,
585 & "backbone fragments."
586 write(iout,*) "fragment, weights:"
588 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
589 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
590 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
593 iset=mod(kolor,nset)+1
596 if(me.eq.king.or..not.out1file)
597 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
600 c------------------------------------------------------------------------------
603 C Read molecular data.
605 implicit real*8 (a-h,o-z)
611 include 'COMMON.IOUNITS'
614 include 'COMMON.INTERACT'
615 include 'COMMON.LOCAL'
616 include 'COMMON.NAMES'
617 include 'COMMON.CHAIN'
618 include 'COMMON.FFIELD'
619 include 'COMMON.SBRIDGE'
620 include 'COMMON.HEADER'
621 include 'COMMON.CONTROL'
622 include 'COMMON.DBASE'
623 include 'COMMON.THREAD'
624 include 'COMMON.CONTACTS'
625 include 'COMMON.TORCNSTR'
626 include 'COMMON.TIME1'
627 include 'COMMON.BOUNDS'
629 include 'COMMON.REMD'
630 include 'COMMON.SETUP'
631 character*4 sequence(maxres)
633 double precision x(maxvar)
634 character*256 pdbfile
635 character*320 weightcard
636 character*80 weightcard_t,ucase
637 dimension itype_pdb(maxres)
638 common /pizda/ itype_pdb
639 logical seq_comp,fail
640 double precision energia(0:n_ene)
647 C Read weights of the subsequent energy terms.
660 if(me.eq.king.or..not.out1file) then
661 write (iout,*) 'Reading ',hremd,' sets of weights for HREMD'
662 write (iout,*) 'Current weights for processor ',
663 & me,' set ',i2set(me)
667 call card_concat(weightcard)
668 call reada(weightcard,'WLONG',wlong,1.0D0)
669 call reada(weightcard,'WSC',wsc,wlong)
670 call reada(weightcard,'WSCP',wscp,wlong)
671 call reada(weightcard,'WELEC',welec,1.0D0)
672 call reada(weightcard,'WVDWPP',wvdwpp,welec)
673 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
674 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
675 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
676 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
677 call reada(weightcard,'WTURN3',wturn3,1.0D0)
678 call reada(weightcard,'WTURN4',wturn4,1.0D0)
679 call reada(weightcard,'WTURN6',wturn6,1.0D0)
680 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
681 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
682 call reada(weightcard,'WBOND',wbond,1.0D0)
683 call reada(weightcard,'WTOR',wtor,1.0D0)
684 call reada(weightcard,'WTORD',wtor_d,1.0D0)
685 call reada(weightcard,'WANG',wang,1.0D0)
686 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
687 call reada(weightcard,'SCAL14',scal14,0.4D0)
688 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
689 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
690 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
691 call reada(weightcard,'TEMP0',temp0,300.0d0)
692 if (index(weightcard,'SOFT').gt.0) ipot=6
693 C 12/1/95 Added weight for the multi-body term WCORR
694 call reada(weightcard,'WCORRH',wcorr,1.0D0)
695 if (wcorr4.gt.0.0d0) wcorr=wcorr4
703 hweights(i,7)=wel_loc
706 hweights(i,10)=wturn6
708 hweights(i,12)=wscloc
710 hweights(i,14)=wtor_d
711 hweights(i,15)=wstrain
712 hweights(i,16)=wvdwpp
714 hweights(i,18)=scal14
715 hweights(i,21)=wsccor
720 weights(i)=hweights(i2set(me),i)
744 call card_concat(weightcard)
745 call reada(weightcard,'WLONG',wlong,1.0D0)
746 call reada(weightcard,'WSC',wsc,wlong)
747 call reada(weightcard,'WSCP',wscp,wlong)
748 call reada(weightcard,'WELEC',welec,1.0D0)
749 call reada(weightcard,'WVDWPP',wvdwpp,welec)
750 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
751 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
752 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
753 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
754 call reada(weightcard,'WTURN3',wturn3,1.0D0)
755 call reada(weightcard,'WTURN4',wturn4,1.0D0)
756 call reada(weightcard,'WTURN6',wturn6,1.0D0)
757 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
758 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
759 call reada(weightcard,'WBOND',wbond,1.0D0)
760 call reada(weightcard,'WTOR',wtor,1.0D0)
761 call reada(weightcard,'WTORD',wtor_d,1.0D0)
762 call reada(weightcard,'WANG',wang,1.0D0)
763 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
764 call reada(weightcard,'WDFAD',wdfa_dist,0.0d0)
765 call reada(weightcard,'WDFAT',wdfa_tor,0.0d0)
766 call reada(weightcard,'WDFAN',wdfa_nei,0.0d0)
767 call reada(weightcard,'WDFAB',wdfa_beta,0.0d0)
768 call reada(weightcard,'SCAL14',scal14,0.4D0)
769 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
770 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
771 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
772 call reada(weightcard,'TEMP0',temp0,300.0d0)
773 if (index(weightcard,'SOFT').gt.0) ipot=6
774 C 12/1/95 Added weight for the multi-body term WCORR
775 call reada(weightcard,'WCORRH',wcorr,1.0D0)
776 if (wcorr4.gt.0.0d0) wcorr=wcorr4
797 weights(25)=wdfa_dist
800 weights(28)=wdfa_beta
802 if(me.eq.king.or..not.out1file)
803 & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
804 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
806 & wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta
807 10 format (/'Energy-term weights (unscaled):'//
808 & 'WSCC= ',f10.6,' (SC-SC)'/
809 & 'WSCP= ',f10.6,' (SC-p)'/
810 & 'WELEC= ',f10.6,' (p-p electr)'/
811 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
812 & 'WBOND= ',f10.6,' (stretching)'/
813 & 'WANG= ',f10.6,' (bending)'/
814 & 'WSCLOC= ',f10.6,' (SC local)'/
815 & 'WTOR= ',f10.6,' (torsional)'/
816 & 'WTORD= ',f10.6,' (double torsional)'/
817 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
818 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
819 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
820 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
821 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
822 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
823 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
824 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
825 & 'WTURN6= ',f10.6,' (turns, 6th order)'/
826 & 'WDFA_D= ',f10.6,' (DFA, distance)' /
827 & 'WDFA_T= ',f10.6,' (DFA, torsional)' /
828 & 'WDFA_N= ',f10.6,' (DFA, number of neighbor)' /
829 & 'WDFA_B= ',f10.6,' (DFA, beta formation)')
830 if(me.eq.king.or..not.out1file)then
831 if (wcorr4.gt.0.0d0) then
832 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
833 & 'between contact pairs of peptide groups'
834 write (iout,'(2(a,f5.3/))')
835 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
836 & 'Range of quenching the correlation terms:',2*delt_corr
837 else if (wcorr.gt.0.0d0) then
838 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
839 & 'between contact pairs of peptide groups'
841 write (iout,'(a,f8.3)')
842 & 'Scaling factor of 1,4 SC-p interactions:',scal14
843 write (iout,'(a,f8.3)')
844 & 'General scaling factor of SC-p interactions:',scalscp
846 r0_corr=cutoff_corr-delt_corr
848 aad(i,1)=scalscp*aad(i,1)
849 aad(i,2)=scalscp*aad(i,2)
850 bad(i,1)=scalscp*bad(i,1)
851 bad(i,2)=scalscp*bad(i,2)
853 call rescale_weights(t_bath)
854 if(me.eq.king.or..not.out1file)
855 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
856 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
858 & wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta
859 22 format (/'Energy-term weights (scaled):'//
860 & 'WSCC= ',f10.6,' (SC-SC)'/
861 & 'WSCP= ',f10.6,' (SC-p)'/
862 & 'WELEC= ',f10.6,' (p-p electr)'/
863 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
864 & 'WBOND= ',f10.6,' (stretching)'/
865 & 'WANG= ',f10.6,' (bending)'/
866 & 'WSCLOC= ',f10.6,' (SC local)'/
867 & 'WTOR= ',f10.6,' (torsional)'/
868 & 'WTORD= ',f10.6,' (double torsional)'/
869 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
870 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
871 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
872 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
873 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
874 & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
875 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
876 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
877 & 'WTURN6= ',f10.6,' (turns, 6th order)'/
878 & 'WDFA_D= ',f10.6,' (DFA, distance)' /
879 & 'WDFA_T= ',f10.6,' (DFA, torsional)' /
880 & 'WDFA_N= ',f10.6,' (DFA, number of neighbor)' /
881 & 'WDFA_B= ',f10.6,' (DFA, beta formation)')
882 if(me.eq.king.or..not.out1file)
883 & write (iout,*) "Reference temperature for weights calculation:",
885 call reada(weightcard,"D0CM",d0cm,3.78d0)
886 call reada(weightcard,"AKCM",akcm,15.1d0)
887 call reada(weightcard,"AKTH",akth,11.0d0)
888 call reada(weightcard,"AKCT",akct,12.0d0)
889 call reada(weightcard,"V1SS",v1ss,-1.08d0)
890 call reada(weightcard,"V2SS",v2ss,7.61d0)
891 call reada(weightcard,"V3SS",v3ss,13.7d0)
892 call reada(weightcard,"EBR",ebr,-5.50D0)
893 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
895 dyn_ss_mask(i)=.false.
899 dyn_ssbond_ij(i,j)=1.0d300
902 call reada(weightcard,"HT",Ht,0.0D0)
904 ss_depth=ebr/wsc-0.25*eps(1,1)
905 Ht=Ht/wsc-0.25*eps(1,1)
906 akcm=akcm*wstrain/wsc
907 akth=akth*wstrain/wsc
908 akct=akct*wstrain/wsc
909 v1ss=v1ss*wstrain/wsc
910 v2ss=v2ss*wstrain/wsc
911 v3ss=v3ss*wstrain/wsc
913 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
916 if(me.eq.king.or..not.out1file) then
917 write (iout,*) "Parameters of the SS-bond potential:"
918 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
920 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
921 write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
922 write (iout,*)" HT",Ht
923 print *,'indpdb=',indpdb,' pdbref=',pdbref
925 if (indpdb.gt.0 .or. pdbref) then
926 read(inp,'(a)') pdbfile
927 if(me.eq.king.or..not.out1file)
928 & write (iout,'(2a)') 'PDB data will be read from file ',
929 & pdbfile(:ilen(pdbfile))
930 open(ipdbin,file=pdbfile,status='old',err=33)
932 33 write (iout,'(a)') 'Error opening PDB file.'
935 c print *,'Begin reading pdb data'
944 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
945 & (crefjlee(j,i+nres),j=1,3)
948 c print *,'Finished reading pdb data'
949 if(me.eq.king.or..not.out1file)
950 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
951 & ' nstart_sup=',nstart_sup
953 itype_pdb(i)=itype(i)
957 nct=nstart_sup+nsup-1
958 call contact(.false.,ncont_ref,icont_ref,co)
961 C Following 2 lines for diagnostics; comment out if not needed
962 write (iout,*) "Before sideadd"
964 if(me.eq.king.or..not.out1file)
965 & write(iout,*)'Adding sidechains'
972 do while (fail.and.nsi.le.maxsi)
973 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
976 if(fail) write(iout,*)'Adding sidechain failed for res ',
977 & i,' after ',nsi,' trials'
980 C 10/03/12 Adam: Recalculate coordinates with new side chain positions
983 C Following 2 lines for diagnostics; comment out if not needed
984 c write (iout,*) "After sideadd"
987 if (indpdb.eq.0) then
988 C Read sequence if not taken from the pdb file.
990 c print *,'nres=',nres
991 if (iscode.gt.0) then
992 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
994 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
996 C Convert sequence to numeric code
998 itype(i)=rescode(i,sequence(i),iscode)
1000 C Assign initial virtual bond lengths
1006 vbld(i+nres)=dsc(itype(i))
1007 vbld_inv(i+nres)=dsc_inv(itype(i))
1008 c write (iout,*) "i",i," itype",itype(i),
1009 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
1013 c print '(20i4)',(itype(i),i=1,nres)
1016 if (itype(i).eq.21 .or. itype(i+1).eq.21) then
1018 if (itype(i).eq.21) then
1022 else if (itype(i+1).ne.20) then
1024 else if (itype(i).ne.20) then
1031 if(me.eq.king.or..not.out1file)then
1032 write (iout,*) "ITEL"
1034 write (iout,*) i,itype(i),itel(i)
1036 print *,'Call Read_Bridge.'
1039 C 8/13/98 Set limits to generating the dihedral angles
1044 read (inp,*) ndih_constr
1045 if (ndih_constr.gt.0) then
1047 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
1048 if(me.eq.king.or..not.out1file)then
1050 & 'There are',ndih_constr,' constraints on phi angles.'
1052 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
1056 phi0(i)=deg2rad*phi0(i)
1057 drange(i)=deg2rad*drange(i)
1059 if(me.eq.king.or..not.out1file)
1060 & write (iout,*) 'FTORS',ftors
1063 phibound(1,ii) = phi0(i)-drange(i)
1064 phibound(2,ii) = phi0(i)+drange(i)
1069 if (me.eq.king) then
1071 write (iout,'(a)') 'Boundaries in phi angle sampling:'
1073 write (iout,'(a3,i5,2f10.1)')
1074 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
1080 cd print *,'NNT=',NNT,' NCT=',NCT
1081 if (itype(1).eq.21) nnt=2
1082 if (itype(nres).eq.21) nct=nct-1
1084 C Bartek:READ init_vars
1085 C Initialize variables!
1086 C Juyong:READ read_info
1087 C READ fragment information!!
1088 C both routines should be in dfa.F file!!
1090 if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
1091 & wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
1093 print*, 'init_dfa_vars finished!'
1095 print*, 'read_dfa_info finished!'
1099 if(me.eq.king.or..not.out1file)
1100 & write (iout,'(a,i3)') 'nsup=',nsup
1102 if (nsup.le.(nct-nnt+1)) then
1103 do i=0,nct-nnt+1-nsup
1104 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
1110 & 'Error - sequences to be superposed do not match.'
1113 do i=0,nsup-(nct-nnt+1)
1114 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
1116 nstart_sup=nstart_sup+i
1122 & 'Error - sequences to be superposed do not match.'
1125 if (nsup.eq.0) nsup=nct-nnt
1126 if (nstart_sup.eq.0) nstart_sup=nnt
1127 if (nstart_seq.eq.0) nstart_seq=nnt
1128 if(me.eq.king.or..not.out1file)
1129 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
1130 & ' nstart_seq=',nstart_seq
1132 c--- Zscore rms -------
1133 if (nz_start.eq.0) nz_start=nnt
1134 if (nz_end.eq.0 .and. nsup.gt.0) then
1136 else if (nz_end.eq.0) then
1139 if(me.eq.king.or..not.out1file)then
1140 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1141 write (iout,*) 'IZ_SC=',iz_sc
1143 c----------------------
1146 if (.not.pdbref) then
1147 call read_angles(inp,*38)
1149 38 write (iout,'(a)') 'Error reading reference structure.'
1151 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1152 stop 'Error reading reference structure'
1156 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
1165 call contact(.true.,ncont_ref,icont_ref,co)
1167 if(me.eq.king.or..not.out1file)
1168 & write (iout,*) 'Contact order:',co
1170 if(me.eq.king.or..not.out1file)
1171 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1174 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1176 if(me.eq.king.or..not.out1file)
1177 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1178 & icont_ref(1,i),' ',
1179 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1183 c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1184 if (constr_dist.gt.0) then
1185 call read_dist_constr
1189 if (constr_homology.gt.0) then
1190 call read_constr_homology
1191 if (indpdb.gt.0 .or. pdbref) then
1194 c(j,i)=crefjlee(j,i)
1195 cref(j,i)=crefjlee(j,i)
1200 write (iout,*) "Array C"
1202 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),
1203 & (c(j,i+nres),j=1,3)
1205 write (iout,*) "Array Cref"
1207 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i),j=1,3),
1208 & (cref(j,i+nres),j=1,3)
1211 call int_from_cart1(.false.)
1212 call sc_loc_geom(.false.)
1214 thetaref(i)=theta(i)
1219 dc(j,i)=c(j,i+1)-c(j,i)
1220 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1225 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1226 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1234 if (nhpb.gt.0) call hpb_partition
1235 c write (iout,*) "After read_dist_constr nhpb",nhpb
1237 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1238 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1239 & modecalc.ne.10) then
1240 C If input structure hasn't been supplied from the PDB file read or generate
1242 if (iranconf.eq.0 .and. .not. extconf) then
1243 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1244 & write (iout,'(a)') 'Initial geometry will be read in.'
1246 read (inp,*) time,potE,uconst,t_bath,
1247 & nss,(ihpb(j),jhpb(j),j=1,nss), nn, (qfrag(i),i=1,nn)
1248 read(inp,'(8f10.5)',end=36,err=36)
1249 & ((c(l,k),l=1,3),k=1,nres),
1250 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1251 call int_from_cart1(.false.)
1254 dc(j,i)=c(j,i+1)-c(j,i)
1255 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1259 if (itype(i).ne.10) then
1261 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1262 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1268 call read_angles(inp,*36)
1271 36 write (iout,'(a)') 'Error reading angle file.'
1273 call mpi_finalize( MPI_COMM_WORLD,IERR )
1275 stop 'Error reading angle file.'
1277 else if (extconf) then
1278 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1279 & write (iout,'(a)') 'Extended chain initial geometry.'
1281 theta(i)=90d0*deg2rad
1284 phi(i)=180d0*deg2rad
1287 alph(i)=110d0*deg2rad
1290 omeg(i)=-120d0*deg2rad
1293 if(me.eq.king.or..not.out1file)
1294 & write (iout,'(a)') 'Random-generated initial geometry.'
1298 if (me.eq.king .or. fg_rank.eq.0 .and. (
1299 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1303 call gen_rand_conf(itmp,*30)
1305 30 write (iout,*) 'Failed to generate random conformation',
1306 & ', itrial=',itrial
1307 write (*,*) 'Processor:',me,
1308 & ' Failed to generate random conformation',
1318 write (iout,'(a,i3,a)') 'Processor:',me,
1319 & ' error in generating random conformation.'
1320 write (*,'(a,i3,a)') 'Processor:',me,
1321 & ' error in generating random conformation.'
1324 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1331 elseif (modecalc.eq.4) then
1332 read (inp,'(a)') intinname
1333 open (intin,file=intinname,status='old',err=333)
1334 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1335 & write (iout,'(a)') 'intinname',intinname
1336 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1338 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1340 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1342 stop 'Error opening angle file.'
1346 C Generate distance constraints, if the PDB structure is to be regularized.
1347 if (nthread.gt.0) then
1348 call read_threadbase
1350 write (iout,*) "READRTNS: Calling setup_var"
1353 if (me.eq.king .or. .not. out1file)
1355 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1356 write (iout,'(/a,i3,a)')
1357 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1358 write (iout,'(20i4)') (iss(i),i=1,ns)
1360 write(iout,*)"Running with dynamic disulfide-bond formation"
1362 write (iout,'(/a/)') 'Pre-formed links are:'
1368 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1369 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1375 if (ns.gt.0.and.dyn_ss) then
1379 forcon(i-nss)=forcon(i)
1386 dyn_ss_mask(iss(i))=.true.
1389 if (i2ndstr.gt.0) call secstrp2dihc
1390 c call geom_to_var(nvar,x)
1391 c call etotal(energia(0))
1392 c call enerprint(energia(0))
1393 c call briefout(0,etot)
1395 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1396 cd write (iout,'(a)') 'Variable list:'
1397 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1399 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1400 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1401 & 'Processor',myrank,': end reading molecular data.'
1405 c--------------------------------------------------------------------------
1406 logical function seq_comp(itypea,itypeb,length)
1408 integer length,itypea(length),itypeb(length)
1411 if (itypea(i).ne.itypeb(i)) then
1419 c-----------------------------------------------------------------------------
1420 subroutine read_bridge
1421 C Read information about disulfide bridges.
1422 implicit real*8 (a-h,o-z)
1423 include 'DIMENSIONS'
1427 include 'COMMON.IOUNITS'
1428 include 'COMMON.GEO'
1429 include 'COMMON.VAR'
1430 include 'COMMON.INTERACT'
1431 include 'COMMON.LOCAL'
1432 include 'COMMON.NAMES'
1433 include 'COMMON.CHAIN'
1434 include 'COMMON.FFIELD'
1435 include 'COMMON.SBRIDGE'
1436 include 'COMMON.HEADER'
1437 include 'COMMON.CONTROL'
1438 include 'COMMON.DBASE'
1439 include 'COMMON.THREAD'
1440 include 'COMMON.TIME1'
1441 include 'COMMON.SETUP'
1442 C Read bridging residues.
1443 read (inp,*) ns,(iss(i),i=1,ns)
1445 if(me.eq.king.or..not.out1file)
1446 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1447 C Check whether the specified bridging residues are cystines.
1449 if (itype(iss(i)).ne.1) then
1450 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1451 & 'Do you REALLY think that the residue ',
1452 & restyp(itype(iss(i))),i,
1453 & ' can form a disulfide bridge?!!!'
1454 write (*,'(2a,i3,a)')
1455 & 'Do you REALLY think that the residue ',
1456 & restyp(itype(iss(i))),i,
1457 & ' can form a disulfide bridge?!!!'
1459 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1464 C Read preformed bridges.
1466 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1468 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1471 C Check if the residues involved in bridges are in the specified list of
1472 C bridging residues.
1475 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1476 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1477 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1478 & ' contains residues present in other pairs.'
1479 write (*,'(a,i3,a)') 'Disulfide pair',i,
1480 & ' contains residues present in other pairs.'
1482 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1488 if (ihpb(i).eq.iss(j)) goto 10
1490 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1493 if (jhpb(i).eq.iss(j)) goto 20
1495 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1501 ihpb(i)=ihpb(i)+nres
1502 jhpb(i)=jhpb(i)+nres
1508 c----------------------------------------------------------------------------
1509 subroutine read_x(kanal,*)
1510 implicit real*8 (a-h,o-z)
1511 include 'DIMENSIONS'
1512 include 'COMMON.GEO'
1513 include 'COMMON.VAR'
1514 include 'COMMON.CHAIN'
1515 include 'COMMON.IOUNITS'
1516 include 'COMMON.CONTROL'
1517 include 'COMMON.LOCAL'
1518 include 'COMMON.INTERACT'
1519 c Read coordinates from input
1521 read(kanal,'(8f10.5)',end=10,err=10)
1522 & ((c(l,k),l=1,3),k=1,nres),
1523 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1526 c(j,2*nres)=c(j,nres)
1528 call int_from_cart1(.false.)
1531 dc(j,i)=c(j,i+1)-c(j,i)
1532 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1536 if (itype(i).ne.10) then
1538 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1539 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1547 c----------------------------------------------------------------------------
1548 subroutine read_threadbase
1549 implicit real*8 (a-h,o-z)
1550 include 'DIMENSIONS'
1551 include 'COMMON.IOUNITS'
1552 include 'COMMON.GEO'
1553 include 'COMMON.VAR'
1554 include 'COMMON.INTERACT'
1555 include 'COMMON.LOCAL'
1556 include 'COMMON.NAMES'
1557 include 'COMMON.CHAIN'
1558 include 'COMMON.FFIELD'
1559 include 'COMMON.SBRIDGE'
1560 include 'COMMON.HEADER'
1561 include 'COMMON.CONTROL'
1562 include 'COMMON.DBASE'
1563 include 'COMMON.THREAD'
1564 include 'COMMON.TIME1'
1565 C Read pattern database for threading.
1566 read (icbase,*) nseq
1568 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1569 & nres_base(2,i),nres_base(3,i)
1570 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1572 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1573 c & nres_base(2,i),nres_base(3,i)
1574 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1578 if (weidis.eq.0.0D0) weidis=0.1D0
1587 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1588 write (iout,'(a,i5)') 'nexcl: ',nexcl
1589 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1592 c------------------------------------------------------------------------------
1593 subroutine setup_var
1594 implicit real*8 (a-h,o-z)
1595 include 'DIMENSIONS'
1596 include 'COMMON.IOUNITS'
1597 include 'COMMON.GEO'
1598 include 'COMMON.VAR'
1599 include 'COMMON.INTERACT'
1600 include 'COMMON.LOCAL'
1601 include 'COMMON.NAMES'
1602 include 'COMMON.CHAIN'
1603 include 'COMMON.FFIELD'
1604 include 'COMMON.SBRIDGE'
1605 include 'COMMON.HEADER'
1606 include 'COMMON.CONTROL'
1607 include 'COMMON.DBASE'
1608 include 'COMMON.THREAD'
1609 include 'COMMON.TIME1'
1610 C Set up variable list.
1616 if (itype(i).ne.10) then
1618 ialph(i,1)=nvar+nside
1622 if (indphi.gt.0) then
1624 else if (indback.gt.0) then
1629 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1632 c----------------------------------------------------------------------------
1633 subroutine gen_dist_constr
1634 C Generate CA distance constraints.
1635 implicit real*8 (a-h,o-z)
1636 include 'DIMENSIONS'
1637 include 'COMMON.IOUNITS'
1638 include 'COMMON.GEO'
1639 include 'COMMON.VAR'
1640 include 'COMMON.INTERACT'
1641 include 'COMMON.LOCAL'
1642 include 'COMMON.NAMES'
1643 include 'COMMON.CHAIN'
1644 include 'COMMON.FFIELD'
1645 include 'COMMON.SBRIDGE'
1646 include 'COMMON.HEADER'
1647 include 'COMMON.CONTROL'
1648 include 'COMMON.DBASE'
1649 include 'COMMON.THREAD'
1650 include 'COMMON.TIME1'
1651 dimension itype_pdb(maxres)
1652 common /pizda/ itype_pdb
1654 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1655 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1656 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1658 do i=nstart_sup,nstart_sup+nsup-1
1659 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1660 cd & ' seq_pdb', restyp(itype_pdb(i))
1661 do j=i+2,nstart_sup+nsup-1
1663 ihpb(nhpb)=i+nstart_seq-nstart_sup
1664 jhpb(nhpb)=j+nstart_seq-nstart_sup
1666 dhpb(nhpb)=dist(i,j)
1669 cd write (iout,'(a)') 'Distance constraints:'
1674 cd if (ii.gt.nres) then
1679 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1680 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1681 cd & dhpb(i),forcon(i)
1685 c----------------------------------------------------------------------------
1687 implicit real*8 (a-h,o-z)
1688 include 'DIMENSIONS'
1689 include 'COMMON.MAP'
1690 include 'COMMON.IOUNITS'
1691 character*3 angid(4) /'THE','PHI','ALP','OME'/
1692 character*80 mapcard,ucase
1694 read (inp,'(a)') mapcard
1695 mapcard=ucase(mapcard)
1696 if (index(mapcard,'PHI').gt.0) then
1698 else if (index(mapcard,'THE').gt.0) then
1700 else if (index(mapcard,'ALP').gt.0) then
1702 else if (index(mapcard,'OME').gt.0) then
1705 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1706 stop 'Error - illegal variable spec in MAP card.'
1708 call readi (mapcard,'RES1',res1(imap),0)
1709 call readi (mapcard,'RES2',res2(imap),0)
1710 if (res1(imap).eq.0) then
1711 res1(imap)=res2(imap)
1712 else if (res2(imap).eq.0) then
1713 res2(imap)=res1(imap)
1715 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1717 & 'Error - illegal definition of variable group in MAP.'
1718 stop 'Error - illegal definition of variable group in MAP.'
1720 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1721 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1722 call readi(mapcard,'NSTEP',nstep(imap),0)
1723 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1725 & 'Illegal boundary and/or step size specification in MAP.'
1726 stop 'Illegal boundary and/or step size specification in MAP.'
1731 c----------------------------------------------------------------------------
1732 csa subroutine csaread
1733 csa implicit real*8 (a-h,o-z)
1734 csa include 'DIMENSIONS'
1735 csa include 'COMMON.IOUNITS'
1736 csa include 'COMMON.GEO'
1737 csa include 'COMMON.CSA'
1738 csa include 'COMMON.BANK'
1739 csa include 'COMMON.CONTROL'
1740 csa character*80 ucase
1741 csa character*620 mcmcard
1742 csa call card_concat(mcmcard)
1744 csa call readi(mcmcard,'NCONF',nconf,50)
1745 csa call readi(mcmcard,'NADD',nadd,0)
1746 csa call readi(mcmcard,'JSTART',jstart,1)
1747 csa call readi(mcmcard,'JEND',jend,1)
1748 csa call readi(mcmcard,'NSTMAX',nstmax,500000)
1749 csa call readi(mcmcard,'N0',n0,1)
1750 csa call readi(mcmcard,'N1',n1,6)
1751 csa call readi(mcmcard,'N2',n2,4)
1752 csa call readi(mcmcard,'N3',n3,0)
1753 csa call readi(mcmcard,'N4',n4,0)
1754 csa call readi(mcmcard,'N5',n5,0)
1755 csa call readi(mcmcard,'N6',n6,10)
1756 csa call readi(mcmcard,'N7',n7,0)
1757 csa call readi(mcmcard,'N8',n8,0)
1758 csa call readi(mcmcard,'N9',n9,0)
1759 csa call readi(mcmcard,'N14',n14,0)
1760 csa call readi(mcmcard,'N15',n15,0)
1761 csa call readi(mcmcard,'N16',n16,0)
1762 csa call readi(mcmcard,'N17',n17,0)
1763 csa call readi(mcmcard,'N18',n18,0)
1765 csa vdisulf=(index(mcmcard,'DYNSS').gt.0)
1767 csa call readi(mcmcard,'NDIFF',ndiff,2)
1768 csa call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1769 csa call readi(mcmcard,'IS1',is1,1)
1770 csa call readi(mcmcard,'IS2',is2,8)
1771 csa call readi(mcmcard,'NRAN0',nran0,4)
1772 csa call readi(mcmcard,'NRAN1',nran1,2)
1773 csa call readi(mcmcard,'IRR',irr,1)
1774 csa call readi(mcmcard,'NSEED',nseed,20)
1775 csa call readi(mcmcard,'NTOTAL',ntotal,10000)
1776 csa call reada(mcmcard,'CUT1',cut1,2.0d0)
1777 csa call reada(mcmcard,'CUT2',cut2,5.0d0)
1778 csa call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1779 csa call readi(mcmcard,'ICMAX',icmax,3)
1780 csa call readi(mcmcard,'IRESTART',irestart,0)
1781 csac!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1784 csa call reada(mcmcard,'DELE',dele,20.0d0)
1785 csa call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1786 csa call readi(mcmcard,'IREF',iref,0)
1787 csa call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1788 csa call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1789 csa call readi(mcmcard,'NCONF_IN',nconf_in,0)
1790 csa call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1791 csa write (iout,*) "NCONF_IN",nconf_in
1794 c----------------------------------------------------------------------------
1795 cfmc subroutine mcmfread
1796 cfmc implicit real*8 (a-h,o-z)
1797 cfmc include 'DIMENSIONS'
1798 cfmc include 'COMMON.MCMF'
1799 cfmc include 'COMMON.IOUNITS'
1800 cfmc include 'COMMON.GEO'
1801 cfmc character*80 ucase
1802 cfmc character*620 mcmcard
1803 cfmc call card_concat(mcmcard)
1805 cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
1806 cfmc write(iout,*)'MAXRANT=',maxrant
1807 cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1808 cfmc write(iout,*)'MAXFAM=',maxfam
1809 cfmc call readi(mcmcard,'NNET1',nnet1,5)
1810 cfmc write(iout,*)'NNET1=',nnet1
1811 cfmc call readi(mcmcard,'NNET2',nnet2,4)
1812 cfmc write(iout,*)'NNET2=',nnet2
1813 cfmc call readi(mcmcard,'NNET3',nnet3,4)
1814 cfmc write(iout,*)'NNET3=',nnet3
1815 cfmc call readi(mcmcard,'ILASTT',ilastt,0)
1816 cfmc write(iout,*)'ILASTT=',ilastt
1817 cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1818 cfmc write(iout,*)'MAXSTR=',maxstr
1819 cfmc maxstr_f=maxstr/maxfam
1820 cfmc write(iout,*)'MAXSTR_F=',maxstr_f
1821 cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
1822 cfmc write(iout,*)'NMCMF=',nmcmf
1823 cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1824 cfmc write(iout,*)'IFOCUS=',ifocus
1825 cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1826 cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
1827 cfmc call readi(mcmcard,'INTPRT',intprt,1000)
1828 cfmc write(iout,*)'INTPRT=',intprt
1829 cfmc call readi(mcmcard,'IPRT',iprt,100)
1830 cfmc write(iout,*)'IPRT=',iprt
1831 cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
1832 cfmc write(iout,*)'IMAXTR=',imaxtr
1833 cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
1834 cfmc write(iout,*)'MAXEVEN=',maxeven
1835 cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1836 cfmc write(iout,*)'MAXEVEN1=',maxeven1
1837 cfmc call readi(mcmcard,'INIMIN',inimin,200)
1838 cfmc write(iout,*)'INIMIN=',inimin
1839 cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1840 cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
1841 cfmc call readi(mcmcard,'NTHREAD',nthread,5)
1842 cfmc write(iout,*)'NTHREAD=',nthread
1843 cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1844 cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1845 cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
1846 cfmc write(iout,*)'MAXPERT=',maxpert
1847 cfmc call readi(mcmcard,'IRMSD',irmsd,1)
1848 cfmc write(iout,*)'IRMSD=',irmsd
1849 cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1850 cfmc write(iout,*)'DENEMIN=',denemin
1851 cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1852 cfmc write(iout,*)'RCUT1S=',rcut1s
1853 cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1854 cfmc write(iout,*)'RCUT1E=',rcut1e
1855 cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1856 cfmc write(iout,*)'RCUT2S=',rcut2s
1857 cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1858 cfmc write(iout,*)'RCUT2E=',rcut2e
1859 cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1860 cfmc write(iout,*)'DPERT1=',d_pert1
1861 cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1862 cfmc write(iout,*)'DPERT1A=',d_pert1a
1863 cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1864 cfmc write(iout,*)'DPERT2=',d_pert2
1865 cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1866 cfmc write(iout,*)'DPERT2A=',d_pert2a
1867 cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1868 cfmc write(iout,*)'DPERT2B=',d_pert2b
1869 cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1870 cfmc write(iout,*)'DPERT2C=',d_pert2c
1871 cfmc d_pert1=deg2rad*d_pert1
1872 cfmc d_pert1a=deg2rad*d_pert1a
1873 cfmc d_pert2=deg2rad*d_pert2
1874 cfmc d_pert2a=deg2rad*d_pert2a
1875 cfmc d_pert2b=deg2rad*d_pert2b
1876 cfmc d_pert2c=deg2rad*d_pert2c
1877 cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1878 cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
1879 cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1880 cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
1881 cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1882 cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1883 cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1884 cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1885 cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1886 cfmc write(iout,*)'RCUTINI=',rcutini
1887 cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
1888 cfmc write(iout,*)'GRAT=',grat
1889 cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1890 cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
1894 c----------------------------------------------------------------------------
1896 implicit real*8 (a-h,o-z)
1897 include 'DIMENSIONS'
1898 include 'COMMON.MCM'
1899 include 'COMMON.MCE'
1900 include 'COMMON.IOUNITS'
1902 character*320 mcmcard
1903 call card_concat(mcmcard)
1904 call readi(mcmcard,'MAXACC',maxacc,100)
1905 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1906 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1907 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1908 call readi(mcmcard,'MAXREPM',maxrepm,200)
1909 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1910 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1911 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1912 call reada(mcmcard,'E_UP',e_up,5.0D0)
1913 call reada(mcmcard,'DELTE',delte,0.1D0)
1914 call readi(mcmcard,'NSWEEP',nsweep,5)
1915 call readi(mcmcard,'NSTEPH',nsteph,0)
1916 call readi(mcmcard,'NSTEPC',nstepc,0)
1917 call reada(mcmcard,'TMIN',tmin,298.0D0)
1918 call reada(mcmcard,'TMAX',tmax,298.0D0)
1919 call readi(mcmcard,'NWINDOW',nwindow,0)
1920 call readi(mcmcard,'PRINT_MC',print_mc,0)
1921 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1922 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1923 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1924 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1925 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1926 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1927 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1928 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1929 if (nwindow.gt.0) then
1930 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1932 winlen(i)=winend(i)-winstart(i)+1
1935 if (tmax.lt.tmin) tmax=tmin
1936 if (tmax.eq.tmin) then
1940 if (nstepc.gt.0 .and. nsteph.gt.0) then
1941 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1942 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1944 C Probabilities of different move types
1945 sumpro_type(0)=0.0D0
1946 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1947 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1948 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1949 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1950 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1951 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1952 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1954 print *,'i',i,' sumprotype',sumpro_type(i)
1955 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1956 print *,'i',i,' sumprotype',sumpro_type(i)
1960 c----------------------------------------------------------------------------
1961 subroutine read_minim
1962 implicit real*8 (a-h,o-z)
1963 include 'DIMENSIONS'
1964 include 'COMMON.MINIM'
1965 include 'COMMON.IOUNITS'
1967 character*320 minimcard
1968 call card_concat(minimcard)
1969 call readi(minimcard,'MAXMIN',maxmin,2000)
1970 call readi(minimcard,'MAXFUN',maxfun,5000)
1971 call readi(minimcard,'MINMIN',minmin,maxmin)
1972 call readi(minimcard,'MINFUN',minfun,maxmin)
1973 call reada(minimcard,'TOLF',tolf,1.0D-2)
1974 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1975 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1976 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1977 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1978 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1979 & 'Options in energy minimization:'
1980 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1981 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1982 & 'MinMin:',MinMin,' MinFun:',MinFun,
1983 & ' TolF:',TolF,' RTolF:',RTolF
1986 c----------------------------------------------------------------------------
1987 subroutine read_angles(kanal,*)
1988 implicit real*8 (a-h,o-z)
1989 include 'DIMENSIONS'
1990 include 'COMMON.GEO'
1991 include 'COMMON.VAR'
1992 include 'COMMON.CHAIN'
1993 include 'COMMON.IOUNITS'
1994 include 'COMMON.CONTROL'
1995 c Read angles from input
1997 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1998 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1999 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
2000 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
2003 c 9/7/01 avoid 180 deg valence angle
2004 if (theta(i).gt.179.99d0) theta(i)=179.99d0
2006 theta(i)=deg2rad*theta(i)
2007 phi(i)=deg2rad*phi(i)
2008 alph(i)=deg2rad*alph(i)
2009 omeg(i)=deg2rad*omeg(i)
2014 c----------------------------------------------------------------------------
2015 subroutine reada(rekord,lancuch,wartosc,default)
2017 character*(*) rekord,lancuch
2018 double precision wartosc,default
2021 iread=index(rekord,lancuch)
2022 if (iread.eq.0) then
2026 iread=iread+ilen(lancuch)+1
2027 read (rekord(iread:),*,err=10,end=10) wartosc
2032 c----------------------------------------------------------------------------
2033 subroutine readi(rekord,lancuch,wartosc,default)
2035 character*(*) rekord,lancuch
2036 integer wartosc,default
2039 iread=index(rekord,lancuch)
2040 if (iread.eq.0) then
2044 iread=iread+ilen(lancuch)+1
2045 read (rekord(iread:),*,err=10,end=10) wartosc
2050 c----------------------------------------------------------------------------
2051 subroutine multreadi(rekord,lancuch,tablica,dim,default)
2054 integer tablica(dim),default
2055 character*(*) rekord,lancuch
2062 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2063 if (iread.eq.0) return
2064 iread=iread+ilen(lancuch)+1
2065 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2068 c----------------------------------------------------------------------------
2069 subroutine multreada(rekord,lancuch,tablica,dim,default)
2072 double precision tablica(dim),default
2073 character*(*) rekord,lancuch
2080 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2081 if (iread.eq.0) return
2082 iread=iread+ilen(lancuch)+1
2083 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2086 c----------------------------------------------------------------------------
2087 subroutine openunits
2088 implicit real*8 (a-h,o-z)
2089 include 'DIMENSIONS'
2092 character*16 form,nodename
2095 include 'COMMON.SETUP'
2096 include 'COMMON.IOUNITS'
2098 include 'COMMON.CONTROL'
2099 integer lenpre,lenpot,ilen,lentmp
2101 character*3 out1file_text,ucase
2104 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
2105 call getenv_loc("PREFIX",prefix)
2107 call getenv_loc("POT",pot)
2108 call getenv_loc("DIRTMP",tmpdir)
2109 call getenv_loc("CURDIR",curdir)
2110 call getenv_loc("OUT1FILE",out1file_text)
2111 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
2112 out1file_text=ucase(out1file_text)
2113 if (out1file_text(1:1).eq."Y") then
2116 out1file=fg_rank.gt.0
2121 if (lentmp.gt.0) then
2122 write (*,'(80(1h!))')
2123 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
2124 write (*,'(80(1h!))')
2125 write (*,*)"All output files will be on node /tmp directory."
2127 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
2128 if (me.eq.king) then
2129 write (*,*) "The master node is ",nodename
2130 else if (fg_rank.eq.0) then
2131 write (*,*) "I am the CG slave node ",nodename
2133 write (*,*) "I am the FG slave node ",nodename
2136 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
2137 lenpre = lentmp+lenpre+1
2139 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
2140 C Get the names and open the input files
2141 #if defined(WINIFL) || defined(WINPGI)
2142 open(1,file=pref_orig(:ilen(pref_orig))//
2143 & '.inp',status='old',readonly,shared)
2144 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2145 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2146 C Get parameter filenames and open the parameter files.
2147 call getenv_loc('BONDPAR',bondname)
2148 open (ibond,file=bondname,status='old',readonly,shared)
2149 call getenv_loc('THETPAR',thetname)
2150 open (ithep,file=thetname,status='old',readonly,shared)
2152 call getenv_loc('THETPARPDB',thetname_pdb)
2153 open (ithep_pdb,file=thetname_pdb,status='old',readonly,shared)
2155 call getenv_loc('ROTPAR',rotname)
2156 open (irotam,file=rotname,status='old',readonly,shared)
2158 call getenv_loc('ROTPARPDB',rotname_pdb)
2159 open (irotam_pdb,file=rotname_pdb,status='old',readonly,shared)
2161 call getenv_loc('TORPAR',torname)
2162 open (itorp,file=torname,status='old',readonly,shared)
2163 call getenv_loc('TORDPAR',tordname)
2164 open (itordp,file=tordname,status='old',readonly,shared)
2165 call getenv_loc('FOURIER',fouriername)
2166 open (ifourier,file=fouriername,status='old',readonly,shared)
2167 call getenv_loc('ELEPAR',elename)
2168 open (ielep,file=elename,status='old',readonly,shared)
2169 call getenv_loc('SIDEPAR',sidename)
2170 open (isidep,file=sidename,status='old',readonly,shared)
2171 #elif (defined CRAY) || (defined AIX)
2172 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2174 c print *,"Processor",myrank," opened file 1"
2175 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2176 c print *,"Processor",myrank," opened file 9"
2177 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2178 C Get parameter filenames and open the parameter files.
2179 call getenv_loc('BONDPAR',bondname)
2180 open (ibond,file=bondname,status='old',action='read')
2181 c print *,"Processor",myrank," opened file IBOND"
2182 call getenv_loc('THETPAR',thetname)
2183 open (ithep,file=thetname,status='old',action='read')
2184 c print *,"Processor",myrank," opened file ITHEP"
2186 call getenv_loc('THETPARPDB',thetname_pdb)
2187 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2189 call getenv_loc('ROTPAR',rotname)
2190 open (irotam,file=rotname,status='old',action='read')
2191 c print *,"Processor",myrank," opened file IROTAM"
2193 call getenv_loc('ROTPARPDB',rotname_pdb)
2194 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2196 call getenv_loc('TORPAR',torname)
2197 open (itorp,file=torname,status='old',action='read')
2198 c print *,"Processor",myrank," opened file ITORP"
2199 call getenv_loc('TORDPAR',tordname)
2200 open (itordp,file=tordname,status='old',action='read')
2201 c print *,"Processor",myrank," opened file ITORDP"
2202 call getenv_loc('SCCORPAR',sccorname)
2203 open (isccor,file=sccorname,status='old',action='read')
2204 c print *,"Processor",myrank," opened file ISCCOR"
2205 call getenv_loc('FOURIER',fouriername)
2206 open (ifourier,file=fouriername,status='old',action='read')
2207 c print *,"Processor",myrank," opened file IFOURIER"
2208 call getenv_loc('ELEPAR',elename)
2209 open (ielep,file=elename,status='old',action='read')
2210 c print *,"Processor",myrank," opened file IELEP"
2211 call getenv_loc('SIDEPAR',sidename)
2212 open (isidep,file=sidename,status='old',action='read')
2213 c print *,"Processor",myrank," opened file ISIDEP"
2214 c print *,"Processor",myrank," opened parameter files"
2216 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2217 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2218 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2219 C Get parameter filenames and open the parameter files.
2220 call getenv_loc('BONDPAR',bondname)
2221 open (ibond,file=bondname,status='old')
2222 call getenv_loc('THETPAR',thetname)
2223 open (ithep,file=thetname,status='old')
2225 call getenv_loc('THETPARPDB',thetname_pdb)
2226 open (ithep_pdb,file=thetname_pdb,status='old')
2228 call getenv_loc('ROTPAR',rotname)
2229 open (irotam,file=rotname,status='old')
2231 call getenv_loc('ROTPARPDB',rotname_pdb)
2232 open (irotam_pdb,file=rotname_pdb,status='old')
2234 call getenv_loc('TORPAR',torname)
2235 open (itorp,file=torname,status='old')
2236 call getenv_loc('TORDPAR',tordname)
2237 open (itordp,file=tordname,status='old')
2238 call getenv_loc('SCCORPAR',sccorname)
2239 open (isccor,file=sccorname,status='old')
2240 call getenv_loc('FOURIER',fouriername)
2241 open (ifourier,file=fouriername,status='old')
2242 call getenv_loc('ELEPAR',elename)
2243 open (ielep,file=elename,status='old')
2244 call getenv_loc('SIDEPAR',sidename)
2245 open (isidep,file=sidename,status='old')
2247 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2249 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2250 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2251 C Get parameter filenames and open the parameter files.
2252 call getenv_loc('BONDPAR',bondname)
2253 open (ibond,file=bondname,status='old',action='read')
2254 call getenv_loc('THETPAR',thetname)
2255 open (ithep,file=thetname,status='old',action='read')
2257 call getenv_loc('THETPARPDB',thetname_pdb)
2258 print *,"thetname_pdb ",thetname_pdb
2259 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2260 print *,ithep_pdb," opened"
2262 call getenv_loc('ROTPAR',rotname)
2263 open (irotam,file=rotname,status='old',action='read')
2265 call getenv_loc('ROTPARPDB',rotname_pdb)
2266 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2268 call getenv_loc('TORPAR',torname)
2269 open (itorp,file=torname,status='old',action='read')
2270 call getenv_loc('TORDPAR',tordname)
2271 open (itordp,file=tordname,status='old',action='read')
2272 call getenv_loc('SCCORPAR',sccorname)
2273 open (isccor,file=sccorname,status='old',action='read')
2274 call getenv_loc('FOURIER',fouriername)
2275 open (ifourier,file=fouriername,status='old',action='read')
2276 call getenv_loc('ELEPAR',elename)
2277 open (ielep,file=elename,status='old',action='read')
2278 call getenv_loc('SIDEPAR',sidename)
2279 open (isidep,file=sidename,status='old',action='read')
2283 C 8/9/01 In the newest version SCp interaction constants are read from a file
2284 C Use -DOLDSCP to use hard-coded constants instead.
2286 call getenv_loc('SCPPAR',scpname)
2287 #if defined(WINIFL) || defined(WINPGI)
2288 open (iscpp,file=scpname,status='old',readonly,shared)
2289 #elif (defined CRAY) || (defined AIX)
2290 open (iscpp,file=scpname,status='old',action='read')
2292 open (iscpp,file=scpname,status='old')
2294 open (iscpp,file=scpname,status='old',action='read')
2297 call getenv_loc('PATTERN',patname)
2298 #if defined(WINIFL) || defined(WINPGI)
2299 open (icbase,file=patname,status='old',readonly,shared)
2300 #elif (defined CRAY) || (defined AIX)
2301 open (icbase,file=patname,status='old',action='read')
2303 open (icbase,file=patname,status='old')
2305 open (icbase,file=patname,status='old',action='read')
2308 C Open output file only for CG processes
2309 c print *,"Processor",myrank," fg_rank",fg_rank
2310 if (fg_rank.eq.0) then
2312 if (nodes.eq.1) then
2315 npos = dlog10(dfloat(nodes-1))+1
2317 if (npos.lt.3) npos=3
2318 write (liczba,'(i1)') npos
2319 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2321 write (liczba,form) me
2322 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2323 & liczba(:ilen(liczba))
2324 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2326 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2328 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2329 & liczba(:ilen(liczba))//'.mol2'
2330 cartname=prefix(:lenpre)//'_'//pot(:lenpot)//
2331 & liczba(:ilen(liczba))//'.x'
2332 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2333 & liczba(:ilen(liczba))//'.stat'
2335 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2336 & //liczba(:ilen(liczba))//'.stat')
2337 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2340 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2341 & liczba(:ilen(liczba))//'.const'
2346 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2347 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2348 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2349 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2350 cartname=prefix(:lenpre)//'_'//pot(:lenpot)//'.x'
2351 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2353 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2355 rest2name=prefix(:ilen(prefix))//'.rst'
2357 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2360 #if defined(AIX) || defined(PGI)
2361 if (me.eq.king .or. .not. out1file)
2362 & open(iout,file=outname,status='unknown')
2365 if (fg_rank.gt.0) then
2366 write (liczba,'(i3.3)') myrank/nfgtasks
2367 write (ll,'(bz,i3.3)') fg_rank
2368 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2374 open(igeom,file=intname,status='unknown',position='append')
2375 open(ipdb,file=pdbname,status='unknown')
2376 open(imol2,file=mol2name,status='unknown')
2377 open(istat,file=statname,status='unknown',position='append')
2379 c1out open(iout,file=outname,status='unknown')
2382 if (me.eq.king .or. .not.out1file)
2383 & open(iout,file=outname,status='unknown')
2386 if (fg_rank.gt.0) then
2387 print "Processor",fg_rank," opening output file"
2388 write (liczba,'(i3.3)') myrank/nfgtasks
2389 write (ll,'(bz,i3.3)') fg_rank
2390 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2396 open(igeom,file=intname,status='unknown',access='append')
2397 open(ipdb,file=pdbname,status='unknown')
2398 open(imol2,file=mol2name,status='unknown')
2399 open(istat,file=statname,status='unknown',access='append')
2401 c1out open(iout,file=outname,status='unknown')
2404 csa csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2405 csa csa_seed=prefix(:lenpre)//'.CSA.seed'
2406 csa csa_history=prefix(:lenpre)//'.CSA.history'
2407 csa csa_bank=prefix(:lenpre)//'.CSA.bank'
2408 csa csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2409 csa csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2410 csa csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2411 csac!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2412 csa csa_int=prefix(:lenpre)//'.int'
2413 csa csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2414 csa csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2415 csa csa_in=prefix(:lenpre)//'.CSA.in'
2416 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2419 write (iout,'(80(1h-))')
2420 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2421 write (iout,'(80(1h-))')
2422 write (iout,*) "Input file : ",
2423 & pref_orig(:ilen(pref_orig))//'.inp'
2424 write (iout,*) "Output file : ",
2425 & outname(:ilen(outname))
2427 write (iout,*) "Sidechain potential file : ",
2428 & sidename(:ilen(sidename))
2430 write (iout,*) "SCp potential file : ",
2431 & scpname(:ilen(scpname))
2433 write (iout,*) "Electrostatic potential file : ",
2434 & elename(:ilen(elename))
2435 write (iout,*) "Cumulant coefficient file : ",
2436 & fouriername(:ilen(fouriername))
2437 write (iout,*) "Torsional parameter file : ",
2438 & torname(:ilen(torname))
2439 write (iout,*) "Double torsional parameter file : ",
2440 & tordname(:ilen(tordname))
2441 write (iout,*) "SCCOR parameter file : ",
2442 & sccorname(:ilen(sccorname))
2443 write (iout,*) "Bond & inertia constant file : ",
2444 & bondname(:ilen(bondname))
2445 write (iout,*) "Bending parameter file : ",
2446 & thetname(:ilen(thetname))
2447 write (iout,*) "Rotamer parameter file : ",
2448 & rotname(:ilen(rotname))
2449 write (iout,*) "Threading database : ",
2450 & patname(:ilen(patname))
2452 &write (iout,*)" DIRTMP : ",
2454 write (iout,'(80(1h-))')
2458 c----------------------------------------------------------------------------
2459 subroutine card_concat(card)
2460 implicit real*8 (a-h,o-z)
2461 include 'DIMENSIONS'
2462 include 'COMMON.IOUNITS'
2464 character*80 karta,ucase
2466 read (inp,'(a)') karta
2469 do while (karta(80:80).eq.'&')
2470 card=card(:ilen(card)+1)//karta(:79)
2471 read (inp,'(a)') karta
2474 card=card(:ilen(card)+1)//karta
2477 c----------------------------------------------------------------------------------
2479 implicit real*8 (a-h,o-z)
2480 include 'DIMENSIONS'
2481 include 'COMMON.CHAIN'
2482 include 'COMMON.IOUNITS'
2484 include 'COMMON.CONTROL'
2485 open(irest2,file=rest2name,status='unknown')
2486 read(irest2,*) totT,EK,potE,totE,t_bath
2488 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2491 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2493 if(usampl.or.homol_nset.gt.1) then
2494 read (irest2,*) iset
2499 c---------------------------------------------------------------------------------
2500 subroutine read_fragments
2501 implicit real*8 (a-h,o-z)
2502 include 'DIMENSIONS'
2506 include 'COMMON.SETUP'
2507 include 'COMMON.CHAIN'
2508 include 'COMMON.IOUNITS'
2510 include 'COMMON.CONTROL'
2512 read(inp,*) nset,nfrag,npair,nfrag_back
2513 if(me.eq.king.or..not.out1file)
2514 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2515 & " nfrag_back",nfrag_back
2517 read(inp,*) mset(iset1)
2519 read(inp,*) wfrag(i,iset1),ifrag(1,i,iset1),ifrag(2,i,iset1),
2521 if(me.eq.king.or..not.out1file)
2522 & write(iout,*) "R ",i,wfrag(i,iset1),ifrag(1,i,iset1),
2523 & ifrag(2,i,iset1), qinfrag(i,iset1)
2526 read(inp,*) wpair(i,iset1),ipair(1,i,iset1),ipair(2,i,iset1),
2528 if(me.eq.king.or..not.out1file)
2529 & write(iout,*) "R ",i,wpair(i,iset1),ipair(1,i,iset1),
2530 & ipair(2,i,iset1), qinpair(i,iset1)
2533 read(inp,*) wfrag_back(1,i,iset1),wfrag_back(2,i,iset1),
2534 & wfrag_back(3,i,iset1),
2535 & ifrag_back(1,i,iset1),ifrag_back(2,i,iset1)
2536 if(me.eq.king.or..not.out1file)
2537 & write(iout,*) "A",i,wfrag_back(1,i,iset1),
2538 & wfrag_back(2,i,iset1),
2539 & wfrag_back(3,i,iset1),ifrag_back(1,i,iset1),
2540 & ifrag_back(2,i,iset1)
2545 c-------------------------------------------------------------------------------
2546 subroutine read_dist_constr
2547 implicit real*8 (a-h,o-z)
2548 include 'DIMENSIONS'
2552 include 'COMMON.SETUP'
2553 include 'COMMON.CONTROL'
2554 include 'COMMON.CHAIN'
2555 include 'COMMON.IOUNITS'
2556 include 'COMMON.SBRIDGE'
2557 integer ifrag_(2,100),ipair_(2,100)
2558 double precision wfrag_(100),wpair_(100)
2559 character*500 controlcard
2560 c write (iout,*) "Calling read_dist_constr"
2561 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2563 call card_concat(controlcard)
2564 call readi(controlcard,"NFRAG",nfrag_,0)
2565 call readi(controlcard,"NPAIR",npair_,0)
2566 call readi(controlcard,"NDIST",ndist_,0)
2567 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2568 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2569 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2570 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2571 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2572 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2573 c write (iout,*) "IFRAG"
2575 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2577 c write (iout,*) "IPAIR"
2579 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2581 if (.not.refstr .and. nfrag.gt.0) then
2583 & "ERROR: no reference structure to compute distance restraints"
2585 & "Restraints must be specified explicitly (NDIST=number)"
2588 if (nfrag.lt.2 .and. npair.gt.0) then
2589 write (iout,*) "ERROR: Less than 2 fragments specified",
2590 & " but distance restraints between pairs requested"
2595 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2596 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2597 & ifrag_(2,i)=nstart_sup+nsup-1
2598 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2600 if (wfrag_(i).gt.0.0d0) then
2601 do j=ifrag_(1,i),ifrag_(2,i)-1
2602 do k=j+1,ifrag_(2,i)
2603 c write (iout,*) "j",j," k",k
2605 if (constr_dist.eq.1) then
2610 forcon(nhpb)=wfrag_(i)
2611 else if (constr_dist.eq.2) then
2612 if (ddjk.le.dist_cut) then
2617 forcon(nhpb)=wfrag_(i)
2624 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2627 if (.not.out1file .or. me.eq.king)
2628 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2629 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2631 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2632 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2639 if (wpair_(i).gt.0.0d0) then
2647 do j=ifrag_(1,ii),ifrag_(2,ii)
2648 do k=ifrag_(1,jj),ifrag_(2,jj)
2652 forcon(nhpb)=wpair_(i)
2653 dhpb(nhpb)=dist(j,k)
2655 if (.not.out1file .or. me.eq.king)
2656 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2657 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2659 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2660 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2667 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2668 & ibecarb(i),forcon(nhpb+1)
2669 if (forcon(nhpb+1).gt.0.0d0) then
2671 if (ibecarb(i).gt.0) then
2672 ihpb(i)=ihpb(i)+nres
2673 jhpb(i)=jhpb(i)+nres
2675 if (dhpb(nhpb).eq.0.0d0)
2676 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2680 if (.not.out1file .or. me.eq.king) then
2683 write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ",
2684 & i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i)
2692 c-------------------------------------------------------------------------------
2694 subroutine read_constr_homology
2696 include 'DIMENSIONS'
2700 include 'COMMON.SETUP'
2701 include 'COMMON.CONTROL'
2702 include 'COMMON.CHAIN'
2703 include 'COMMON.IOUNITS'
2705 include 'COMMON.GEO'
2706 include 'COMMON.INTERACT'
2708 c For new homol impl
2710 include 'COMMON.VAR'
2713 c double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
2715 c common /przechowalnia/ odl_temp(maxres,maxres,max_template),
2716 c & sigma_odl_temp(maxres,maxres,max_template)
2718 character*24 model_ki_dist, model_ki_angle
2719 character*500 controlcard
2720 integer ki, i, j, k, l
2721 logical lprn /.true./
2723 c FP - Nov. 2014 Temporary specifications for new vars
2725 double precision rescore_tmp,x12,y12,z12
2726 double precision, dimension (max_template,maxres) :: rescore
2727 character*24 pdbfile,tpl_k_rescore
2728 c -----------------------------------------------------------------
2729 c Reading multiple PDB ref structures and calculation of retraints
2730 c not using pre-computed ones stored in files model_ki_{dist,angle}
2732 c -----------------------------------------------------------------
2735 c Alternative: reading from input
2736 call card_concat(controlcard)
2737 call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
2738 call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
2739 call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
2740 call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
2741 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
2743 call readi(controlcard,"HOMOL_NSET",homol_nset,1)
2744 if (homol_nset.gt.1)then
2745 call card_concat(controlcard)
2746 read(controlcard,*) (waga_homology(i),i=1,homol_nset)
2747 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2748 write(iout,*) "iset homology_weight "
2750 write(iout,*) i,waga_homology(i)
2753 iset=mod(kolor,homol_nset)+1
2756 waga_homology(1)=1.0
2759 cd write (iout,*) "nnt",nnt," nct",nct
2771 c Reading HM global scores (prob not required)
2773 c open (4,file="HMscore")
2774 c do k=1,constr_homology
2775 c read (4,*,end=521) hmscore_tmp
2776 c hmscore(k)=hmscore_tmp ! Another transformation can be used
2777 c write(*,*) "Model", k, ":", hmscore(k)
2781 c write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
2782 do k=1,constr_homology
2784 read(inp,'(a)') pdbfile
2785 c Next stament causes error upon compilation (?)
2786 c if(me.eq.king.or. .not. out1file)
2787 c write (iout,'(2a)') 'PDB data will be read from file ',
2788 c & pdbfile(:ilen(pdbfile))
2789 open(ipdbin,file=pdbfile,status='old',err=33)
2791 33 write (iout,'(a)') 'Error opening PDB file.'
2794 c print *,'Begin reading pdb data'
2796 c Files containing res sim or local scores (former containing sigmas)
2799 write(kic2,'(bz,i2.2)') k
2801 tpl_k_rescore="template"//kic2//".sco"
2802 c tpl_k_sigma_odl="template"//kic2//".sigma_odl"
2803 c tpl_k_sigma_dih="template"//kic2//".sigma_dih"
2804 c tpl_k_sigma_theta="template"//kic2//".sigma_theta"
2805 c tpl_k_sigma_d="template"//kic2//".sigma_d"
2810 c Distance restraints
2813 C Copy the coordinates from reference coordinates (?)
2817 c write (iout,*) "c(",j,i,") =",c(j,i)
2821 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
2823 c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
2824 open (ientin,file=tpl_k_rescore,status='old')
2825 do irec=1,maxdim ! loop for reading res sim
2827 rescore(k,irec)=0.0d0
2830 read (ientin,*,end=1401) rescore_tmp
2831 c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
2832 rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
2833 c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
2838 c open (ientin,file=tpl_k_sigma_odl,status='old')
2839 c do irec=1,maxdim ! loop for reading sigma_odl
2840 c read (ientin,*,end=1401) i, j,
2841 c & sigma_odl_temp(i+nnt-1,j+nnt-1,k) ! new variable (?)
2842 c sigma_odl_temp(j+nnt-1,i+nnt-1,k)= ! which purpose?
2843 c & sigma_odl_temp(i+nnt-1,j+nnt-1,k)
2847 if (waga_dist.ne.0.0d0) then
2849 do i = nnt,nct-2 ! right? without parallel.
2850 do j=i+2,nct ! right?
2851 c do i = 1,nres ! alternative for bounds as used to set initial values in orig. read_constr_homology
2852 c do j=i+2,nres ! ibid
2853 c do i = nnt,nct-2 ! alternative for bounds as used to assign dist restraints in orig. read_constr_homology (s. above)
2854 c do j=i+2,nct ! ibid
2856 c write (iout,*) "k",k
2857 c write (iout,*) "i",i," j",j," constr_homology",
2862 c Attempt to replace dist(i,j) by its definition in ...
2867 distal=dsqrt(x12*x12+y12*y12+z12*z12)
2870 c odl(k,ii)=dist(i,j)
2871 c write (iout,*) "dist(",i,j,") =",dist(i,j)
2872 c write (iout,*) "distal = ",distal
2873 c write (iout,*) "odl(",k,ii,") =",odl(k,ii)
2874 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
2875 c & "rescore(",k,j,") =",rescore(k,j)
2877 c Calculation of sigma from res sim
2879 c if (odl(k,ii).le.6.0d0) then
2880 c sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)
2881 c Other functional forms possible depending on odl(k,ii), eg.
2883 if (odl(k,ii).le.dist_cut) then
2884 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j) ! other exprs possible
2885 c sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)
2887 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* ! sigma ~ rescore ~ error
2888 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
2890 c Following expr replaced by a positive exp argument
2891 c sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
2892 c & dexp(-0.5d0*(odl(k,ii)/dist_cut)**2)
2894 c sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)*
2895 c & dexp(-0.5d0*(odl(k,ii)/dist_cut)**2)
2898 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii)) ! rescore ~ error
2899 c sigma_odl(k,ii)=sigma_odl(k,ii)*sigma_odl(k,ii)
2901 c sigma_odl(k,ii)=sigma_odl_temp(i,j,k)* ! new var read from file (?)
2902 c & sigma_odl_temp(i,j,k) ! not inverse because of use of res. similarity
2904 c read (ientin,*) sigma_odl(k,ii) ! 1st variant
2907 c if (constr_homology.gt.0) call homology_partition
2910 c Theta, dihedral and SC retraints
2912 if (waga_angle.gt.0.0d0) then
2913 c open (ientin,file=tpl_k_sigma_dih,status='old')
2914 c do irec=1,maxres-3 ! loop for reading sigma_dih
2915 c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
2916 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
2917 c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
2918 c & sigma_dih(k,i+nnt-1)
2922 do i = nnt+3,nct ! right? without parallel.
2923 c do i=1,nres ! alternative for bounds acc to readpdb?
2924 c do i=1,nres-3 ! alternative for bounds as used to set initial values in orig. read_constr_homology
2925 c do i=idihconstr_start_homo,idihconstr_end_homo ! with FG parallel.
2926 dih(k,i)=phiref(i) ! right?
2927 c read (ientin,*) sigma_dih(k,i) ! original variant
2928 c write (iout,*) "dih(",k,i,") =",dih(k,i)
2929 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
2930 c & "rescore(",k,i-1,") =",rescore(k,i-1),
2931 c & "rescore(",k,i-2,") =",rescore(k,i-2),
2932 c & "rescore(",k,i-3,") =",rescore(k,i-3)
2934 sigma_dih(k,i)=rescore(k,i)+rescore(k,i-1)+
2935 & rescore(k,i-2)+rescore(k,i-3) ! right expression ?
2937 c write (iout,*) "Raw sigmas for dihedral angle restraints"
2938 c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
2939 c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
2940 c rescore(k,i-2)*rescore(k,i-3) ! right expression ?
2941 c Instead of res sim other local measure of b/b str reliability possible
2942 sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
2943 c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
2944 if (i-nnt-2.gt.lim_dih) lim_dih=i-nnt-2 ! right?
2945 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! original when readin i from file
2949 if (waga_theta.gt.0.0d0) then
2950 c open (ientin,file=tpl_k_sigma_theta,status='old')
2951 c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
2952 c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
2953 c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
2954 c & sigma_theta(k,i+nnt-1)
2959 do i = nnt+2,nct ! right? without parallel.
2960 c do i = i=1,nres ! alternative for bounds acc to readpdb?
2961 c do i=ithet_start,ithet_end ! with FG parallel.
2962 thetatpl(k,i)=thetaref(i)
2963 c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
2964 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
2965 c & "rescore(",k,i-1,") =",rescore(k,i-1),
2966 c & "rescore(",k,i-2,") =",rescore(k,i-2)
2967 c read (ientin,*) sigma_theta(k,i) ! 1st variant
2968 sigma_theta(k,i)=rescore(k,i)+rescore(k,i-1)+
2969 & rescore(k,i-2) ! right expression ?
2970 sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
2972 c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
2973 c rescore(k,i-2) ! right expression ?
2974 c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
2975 if (i-nnt-1.gt.lim_theta) lim_theta=i-nnt-1 ! right?
2979 if (waga_d.gt.0.0d0) then
2980 c open (ientin,file=tpl_k_sigma_d,status='old')
2981 c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
2982 c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
2983 c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
2984 c & sigma_d(k,i+nnt-1)
2989 do i = nnt,nct ! right? without parallel.
2990 c do i=2,nres-1 ! alternative for bounds acc to readpdb?
2991 c do i=loc_start,loc_end ! with FG parallel.
2992 if (itype(i).eq.10) goto 1 ! right?
2996 c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
2997 c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
2998 c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
2999 c write(iout,*) "rescore(",k,i,") =",rescore(k,i)
3000 sigma_d(k,i)=rescore(k,i) ! right expression ?
3001 sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
3003 c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
3004 c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
3005 c read (ientin,*) sigma_d(k,i) ! 1st variant
3006 if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1 ! right?
3012 if (waga_dist.ne.0.0d0) lim_odl=ii
3013 if (constr_homology.gt.0) call homology_partition
3014 if (constr_homology.gt.0) call init_int_table
3015 cd write (iout,*) "homology_partition: lim_theta= ",lim_theta,
3016 cd & "lim_xx=",lim_xx
3017 c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
3018 c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
3022 if (.not.lprn) return
3023 cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
3024 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3025 write (iout,*) "Distance restraints from templates"
3027 write(iout,'(3i5,10(2f16.2,4x))') ii,ires_homo(ii),jres_homo(ii),
3028 & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),ki=1,constr_homology)
3030 write (iout,*) "Dihedral angle restraints from templates"
3032 write (iout,'(i5,10(2f8.2,4x))') i,(rad2deg*dih(ki,i),
3033 & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
3035 write (iout,*) "Virtual-bond angle restraints from templates"
3036 do i=nnt+2,lim_theta
3037 write (iout,'(i5,10(2f8.2,4x))') i,(rad2deg*thetatpl(ki,i),
3038 & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
3040 write (iout,*) "SC restraints from templates"
3042 write(iout,'(i5,10(4f8.2,4x))') i,
3043 & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
3044 & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
3047 c -----------------------------------------------------------------
3050 c----------------------------------------------------------------------
3053 subroutine flush(iu)
3058 subroutine flush(iu)
3064 c------------------------------------------------------------------------------
3065 subroutine copy_to_tmp(source)
3066 include "DIMENSIONS"
3067 include "COMMON.IOUNITS"
3068 character*(*) source
3069 character* 256 tmpfile
3073 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
3074 inquire(file=tmpfile,exist=ex)
3076 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
3077 & " to temporary directory..."
3078 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
3079 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
3083 c------------------------------------------------------------------------------
3084 subroutine move_from_tmp(source)
3085 include "DIMENSIONS"
3086 include "COMMON.IOUNITS"
3087 character*(*) source
3090 write (*,*) "Moving ",source(:ilen(source)),
3091 & " from temporary directory to working directory"
3092 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
3093 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
3096 c------------------------------------------------------------------------------
3097 subroutine random_init(seed)
3099 C Initialize random number generator
3101 implicit real*8 (a-h,o-z)
3102 include 'DIMENSIONS'
3108 logical OKRandom, prng_restart
3110 integer iseed_array(4)
3112 include 'COMMON.IOUNITS'
3113 include 'COMMON.TIME1'
3114 include 'COMMON.THREAD'
3115 include 'COMMON.SBRIDGE'
3116 include 'COMMON.CONTROL'
3117 include 'COMMON.MCM'
3118 include 'COMMON.MAP'
3119 include 'COMMON.HEADER'
3120 csa include 'COMMON.CSA'
3121 include 'COMMON.CHAIN'
3122 include 'COMMON.MUCA'
3124 include 'COMMON.FFIELD'
3125 include 'COMMON.SETUP'
3126 iseed=-dint(dabs(seed))
3127 if (iseed.eq.0) then
3128 write (iout,'(/80(1h*)/20x,a/80(1h*))')
3129 & 'Random seed undefined. The program will stop.'
3130 write (*,'(/80(1h*)/20x,a/80(1h*))')
3131 & 'Random seed undefined. The program will stop.'
3133 call mpi_finalize(mpi_comm_world,ierr)
3135 stop 'Bad random seed.'
3138 if (fg_rank.eq.0) then
3142 if(me.eq.king .or. .not. out1file)
3143 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8
3144 write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8
3145 OKRandom = prng_restart(me,iseedi8)
3148 tmp=65536.0d0**(4-i)
3149 iseed_array(i) = dint(seed/tmp)
3150 seed=seed-iseed_array(i)*tmp
3152 if(me.eq.king .or. .not. out1file)
3153 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
3154 & (iseed_array(i),i=1,4)
3155 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
3156 & (iseed_array(i),i=1,4)
3157 OKRandom = prng_restart(me,iseed_array)
3160 r1=ran_number(0.0D0,1.0D0)
3161 if(me.eq.king .or. .not. out1file)
3162 & write (iout,*) 'ran_num',r1
3163 if (r1.lt.0.0d0) OKRandom=.false.
3165 if (.not.OKRandom) then
3166 write (iout,*) 'PRNG IS NOT WORKING!!!'
3167 print *,'PRNG IS NOT WORKING!!!'
3170 call mpi_abort(mpi_comm_world,error_msg,ierr)
3173 write (iout,*) 'too many processors for parallel prng'
3174 write (*,*) 'too many processors for parallel prng'
3182 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)