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 if (modecalc.eq.8) then
32 inquire (file="fort.40",exist=file_exist)
33 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,f10.5)') i-nss,ihpb(i),jhpb(i),forcon(i)
52 c print *,"Processor",myrank," leaves READRTNS"
55 C-------------------------------------------------------------------------------
56 subroutine read_control
60 implicit real*8 (a-h,o-z)
64 logical OKRandom, prng_restart
67 include 'COMMON.IOUNITS'
68 include 'COMMON.TIME1'
69 include 'COMMON.THREAD'
70 include 'COMMON.SBRIDGE'
71 include 'COMMON.CONTROL'
74 include 'COMMON.HEADER'
76 include 'COMMON.CHAIN'
79 include 'COMMON.FFIELD'
80 include 'COMMON.INTERACT'
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,'SYM',symetr,1)
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
181 crc else if (index(controlcard,'ZSCORE').gt.0) then
183 crc ZSCORE is rm from UNRES, modecalc=9 is available
186 cfcm else if (index(controlcard,'MCMF').gt.0) then
188 else if (index(controlcard,'SOFTREG').gt.0) then
190 else if (index(controlcard,'CHECK_BOND').gt.0) then
192 else if (index(controlcard,'TEST').gt.0) then
194 else if (index(controlcard,'MD').gt.0) then
196 else if (index(controlcard,'RE ').gt.0) then
200 lmuca=index(controlcard,'MUCA').gt.0
201 call readi(controlcard,'MUCADYN',mucadyn,0)
202 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
203 if (lmuca .and. (me.eq.king .or. .not.out1file ))
205 write (iout,*) 'MUCADYN=',mucadyn
206 write (iout,*) 'MUCASMOOTH=',muca_smooth
209 iscode=index(controlcard,'ONE_LETTER')
210 indphi=index(controlcard,'PHI')
211 indback=index(controlcard,'BACK')
212 iranconf=index(controlcard,'RAND_CONF')
213 i2ndstr=index(controlcard,'USE_SEC_PRED')
214 gradout=index(controlcard,'GRADOUT').gt.0
215 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
216 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
217 if (me.eq.king .or. .not.out1file )
218 & write (iout,*) "DISTCHAINMAX",distchainmax
220 if(me.eq.king.or..not.out1file)
221 & write (iout,'(2a)') diagmeth(kdiag),
222 & ' routine used to diagonalize matrices.'
225 c--------------------------------------------------------------------------
226 subroutine read_REMDpar
230 implicit real*8 (a-h,o-z)
232 include 'COMMON.IOUNITS'
233 include 'COMMON.TIME1'
236 include 'COMMON.LANGEVIN'
238 include 'COMMON.LANGEVIN.lang0'
240 include 'COMMON.INTERACT'
241 include 'COMMON.NAMES'
243 include 'COMMON.REMD'
244 include 'COMMON.CONTROL'
245 include 'COMMON.SETUP'
247 character*320 controlcard
248 character*3200 controlcard1
249 integer iremd_m_total
251 if(me.eq.king.or..not.out1file)
252 & write (iout,*) "REMD setup"
254 call card_concat(controlcard)
255 call readi(controlcard,"NREP",nrep,3)
256 call readi(controlcard,"NSTEX",nstex,1000)
257 call reada(controlcard,"RETMIN",retmin,10.0d0)
258 call reada(controlcard,"RETMAX",retmax,1000.0d0)
259 mremdsync=(index(controlcard,'SYNC').gt.0)
260 call readi(controlcard,"NSYN",i_sync_step,100)
261 restart1file=(index(controlcard,'REST1FILE').gt.0)
262 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
263 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
264 if(max_cache_traj_use.gt.max_cache_traj)
265 & max_cache_traj_use=max_cache_traj
266 if(me.eq.king.or..not.out1file) then
267 cd if (traj1file) then
268 crc caching is in testing - NTWX is not ignored
269 cd write (iout,*) "NTWX value is ignored"
270 cd write (iout,*) " trajectory is stored to one file by master"
271 cd write (iout,*) " before exchange at NSTEX intervals"
273 write (iout,*) "NREP= ",nrep
274 write (iout,*) "NSTEX= ",nstex
275 write (iout,*) "SYNC= ",mremdsync
276 write (iout,*) "NSYN= ",i_sync_step
277 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
280 if (index(controlcard,'TLIST').gt.0) then
282 call card_concat(controlcard1)
283 read(controlcard1,*) (remd_t(i),i=1,nrep)
284 if(me.eq.king.or..not.out1file)
285 & write (iout,*)'tlist',(remd_t(i),i=1,nrep)
288 if (index(controlcard,'MLIST').gt.0) then
290 call card_concat(controlcard1)
291 read(controlcard1,*) (remd_m(i),i=1,nrep)
292 if(me.eq.king.or..not.out1file) then
293 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
296 iremd_m_total=iremd_m_total+remd_m(i)
298 write (iout,*) 'Total number of replicas ',iremd_m_total
301 if(me.eq.king.or..not.out1file)
302 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
305 c--------------------------------------------------------------------------
306 subroutine read_MDpar
310 implicit real*8 (a-h,o-z)
312 include 'COMMON.IOUNITS'
313 include 'COMMON.TIME1'
316 include 'COMMON.LANGEVIN'
318 include 'COMMON.LANGEVIN.lang0'
320 include 'COMMON.INTERACT'
321 include 'COMMON.NAMES'
323 include 'COMMON.SETUP'
324 include 'COMMON.CONTROL'
325 include 'COMMON.SPLITELE'
327 character*320 controlcard
329 call card_concat(controlcard)
330 call readi(controlcard,"NSTEP",n_timestep,1000000)
331 call readi(controlcard,"NTWE",ntwe,100)
332 call readi(controlcard,"NTWX",ntwx,1000)
333 call reada(controlcard,"DT",d_time,1.0d-1)
334 call reada(controlcard,"DVMAX",dvmax,2.0d1)
335 call reada(controlcard,"DAMAX",damax,1.0d1)
336 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
337 call readi(controlcard,"LANG",lang,0)
338 RESPA = index(controlcard,"RESPA") .gt. 0
339 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
340 ntime_split0=ntime_split
341 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
342 ntime_split0=ntime_split
343 call reada(controlcard,"R_CUT",r_cut,2.0d0)
344 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
345 rest = index(controlcard,"REST").gt.0
346 tbf = index(controlcard,"TBF").gt.0
347 usampl = index(controlcard,"USAMPL").gt.0
349 mdpdb = index(controlcard,"MDPDB").gt.0
350 call reada(controlcard,"T_BATH",t_bath,300.0d0)
351 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
352 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
353 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
354 if (count_reset_moment.eq.0) count_reset_moment=1000000000
355 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
356 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
357 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
358 if (count_reset_vel.eq.0) count_reset_vel=1000000000
359 large = index(controlcard,"LARGE").gt.0
360 print_compon = index(controlcard,"PRINT_COMPON").gt.0
361 rattle = index(controlcard,"RATTLE").gt.0
362 c if performing umbrella sampling, fragments constrained are read from the fragment file
368 if(me.eq.king.or..not.out1file) then
370 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
372 write (iout,'(a)') "The units are:"
373 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
374 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
375 & " acceleration: angstrom/(48.9 fs)**2"
376 write (iout,'(a)') "energy: kcal/mol, temperature: K"
378 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
379 write (iout,'(a60,f10.5,a)')
380 & "Initial time step of numerical integration:",d_time,
382 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
384 write (iout,'(2a,i4,a)')
385 & "A-MTS algorithm used; initial time step for fast-varying",
386 & " short-range forces split into",ntime_split," steps."
387 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
388 & r_cut," lambda",rlamb
390 write (iout,'(2a,f10.5)')
391 & "Maximum acceleration threshold to reduce the time step",
392 & "/increase split number:",damax
393 write (iout,'(2a,f10.5)')
394 & "Maximum predicted energy drift to reduce the timestep",
395 & "/increase split number:",edriftmax
396 write (iout,'(a60,f10.5)')
397 & "Maximum velocity threshold to reduce velocities:",dvmax
398 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
399 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
400 if (rattle) write (iout,'(a60)')
401 & "Rattle algorithm used to constrain the virtual bonds"
405 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
406 call reada(controlcard,"RWAT",rwat,1.4d0)
407 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
408 surfarea=index(controlcard,"SURFAREA").gt.0
409 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
410 if(me.eq.king.or..not.out1file)then
411 write (iout,'(/a,$)') "Langevin dynamics calculation"
414 & " with direct integration of Langevin equations"
415 else if (lang.eq.2) then
416 write (iout,'(a/)') " with TINKER stochasic MD integrator"
417 else if (lang.eq.3) then
418 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
419 else if (lang.eq.4) then
420 write (iout,'(a/)') " in overdamped mode"
422 write (iout,'(//a,i5)')
423 & "=========== ERROR: Unknown Langevin dynamics mode:",lang
426 write (iout,'(a60,f10.5)') "Temperature:",t_bath
427 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
428 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
429 write (iout,'(a60,f10.5)')
430 & "Scaling factor of the friction forces:",scal_fric
431 if (surfarea) write (iout,'(2a,i10,a)')
432 & "Friction coefficients will be scaled by solvent-accessible",
433 & " surface area every",reset_fricmat," steps."
435 c Calculate friction coefficients and bounds of stochastic forces
436 eta=6*pi*cPoise*etawat
437 if(me.eq.king.or..not.out1file)
438 & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
440 gamp=scal_fric*(pstok+rwat)*eta
441 stdfp=dsqrt(2*Rb*t_bath/d_time)
443 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
444 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
446 if(me.eq.king.or..not.out1file)then
447 write (iout,'(/2a/)')
448 & "Radii of site types and friction coefficients and std's of",
449 & " stochastic forces of fully exposed sites"
450 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
452 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
453 & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
457 if(me.eq.king.or..not.out1file)then
458 write (iout,'(a)') "Berendsen bath calculation"
459 write (iout,'(a60,f10.5)') "Temperature:",t_bath
460 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
462 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
463 & count_reset_moment," steps"
465 & write (iout,'(a,i10,a)')
466 & "Velocities will be reset at random every",count_reset_vel,
470 if(me.eq.king.or..not.out1file)
471 & write (iout,'(a31)') "Microcanonical mode calculation"
473 if(me.eq.king.or..not.out1file)then
474 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
476 write(iout,*) "MD running with constraints."
477 write(iout,*) "Equilibration time ", eq_time, " mtus."
478 write(iout,*) "Constraining ", nfrag," fragments."
479 write(iout,*) "Length of each fragment, weight and q0:"
481 write (iout,*) "Set of restraints #",iset
483 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
484 & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
486 write(iout,*) "constraints between ", npair, "fragments."
487 write(iout,*) "constraint pairs, weights and q0:"
489 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
490 & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
492 write(iout,*) "angle constraints within ", nfrag_back,
493 & "backbone fragments."
494 write(iout,*) "fragment, weights:"
496 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
497 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
498 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
501 iset=mod(kolor,nset)+1
504 if(me.eq.king.or..not.out1file)
505 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
508 c------------------------------------------------------------------------------
511 C Read molecular data.
513 implicit real*8 (a-h,o-z)
519 include 'COMMON.IOUNITS'
522 include 'COMMON.INTERACT'
523 include 'COMMON.LOCAL'
524 include 'COMMON.NAMES'
525 include 'COMMON.CHAIN'
526 include 'COMMON.FFIELD'
527 include 'COMMON.SBRIDGE'
528 include 'COMMON.HEADER'
529 include 'COMMON.CONTROL'
530 include 'COMMON.DBASE'
531 include 'COMMON.THREAD'
532 include 'COMMON.CONTACTS'
533 include 'COMMON.TORCNSTR'
534 include 'COMMON.TIME1'
535 include 'COMMON.BOUNDS'
537 include 'COMMON.SETUP'
538 character*4 sequence(maxres)
540 double precision x(maxvar)
541 character*256 pdbfile
542 character*320 weightcard
543 character*80 weightcard_t,ucase
544 dimension itype_pdb(maxres)
545 common /pizda/ itype_pdb
546 logical seq_comp,fail
547 double precision energia(0:n_ene)
553 C Read weights of the subsequent energy terms.
554 call card_concat(weightcard)
555 call reada(weightcard,'WLONG',wlong,1.0D0)
556 call reada(weightcard,'WSC',wsc,wlong)
557 call reada(weightcard,'WSCP',wscp,wlong)
558 call reada(weightcard,'WELEC',welec,1.0D0)
559 call reada(weightcard,'WVDWPP',wvdwpp,welec)
560 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
561 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
562 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
563 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
564 call reada(weightcard,'WTURN3',wturn3,1.0D0)
565 call reada(weightcard,'WTURN4',wturn4,1.0D0)
566 call reada(weightcard,'WTURN6',wturn6,1.0D0)
567 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
568 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
569 call reada(weightcard,'WBOND',wbond,1.0D0)
570 call reada(weightcard,'WTOR',wtor,1.0D0)
571 call reada(weightcard,'WTORD',wtor_d,1.0D0)
572 call reada(weightcard,'WANG',wang,1.0D0)
573 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
574 call reada(weightcard,'SCAL14',scal14,0.4D0)
575 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
576 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
577 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
578 call reada(weightcard,'TEMP0',temp0,300.0d0)
579 if (index(weightcard,'SOFT').gt.0) ipot=6
580 C 12/1/95 Added weight for the multi-body term WCORR
581 call reada(weightcard,'WCORRH',wcorr,1.0D0)
582 if (wcorr4.gt.0.0d0) wcorr=wcorr4
602 if(me.eq.king.or..not.out1file)
603 & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
604 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
606 10 format (/'Energy-term weights (unscaled):'//
607 & 'WSCC= ',f10.6,' (SC-SC)'/
608 & 'WSCP= ',f10.6,' (SC-p)'/
609 & 'WELEC= ',f10.6,' (p-p electr)'/
610 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
611 & 'WBOND= ',f10.6,' (stretching)'/
612 & 'WANG= ',f10.6,' (bending)'/
613 & 'WSCLOC= ',f10.6,' (SC local)'/
614 & 'WTOR= ',f10.6,' (torsional)'/
615 & 'WTORD= ',f10.6,' (double torsional)'/
616 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
617 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
618 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
619 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
620 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
621 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
622 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
623 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
624 & 'WTURN6= ',f10.6,' (turns, 6th order)')
625 if(me.eq.king.or..not.out1file)then
626 if (wcorr4.gt.0.0d0) then
627 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
628 & 'between contact pairs of peptide groups'
629 write (iout,'(2(a,f5.3/))')
630 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
631 & 'Range of quenching the correlation terms:',2*delt_corr
632 else if (wcorr.gt.0.0d0) then
633 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
634 & 'between contact pairs of peptide groups'
636 write (iout,'(a,f8.3)')
637 & 'Scaling factor of 1,4 SC-p interactions:',scal14
638 write (iout,'(a,f8.3)')
639 & 'General scaling factor of SC-p interactions:',scalscp
641 r0_corr=cutoff_corr-delt_corr
643 aad(i,1)=scalscp*aad(i,1)
644 aad(i,2)=scalscp*aad(i,2)
645 bad(i,1)=scalscp*bad(i,1)
646 bad(i,2)=scalscp*bad(i,2)
648 call rescale_weights(t_bath)
649 if(me.eq.king.or..not.out1file)
650 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
651 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
653 22 format (/'Energy-term weights (scaled):'//
654 & 'WSCC= ',f10.6,' (SC-SC)'/
655 & 'WSCP= ',f10.6,' (SC-p)'/
656 & 'WELEC= ',f10.6,' (p-p electr)'/
657 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
658 & 'WBOND= ',f10.6,' (stretching)'/
659 & 'WANG= ',f10.6,' (bending)'/
660 & 'WSCLOC= ',f10.6,' (SC local)'/
661 & 'WTOR= ',f10.6,' (torsional)'/
662 & 'WTORD= ',f10.6,' (double torsional)'/
663 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
664 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
665 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
666 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
667 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
668 & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
669 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
670 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
671 & 'WTURN6= ',f10.6,' (turns, 6th order)')
672 if(me.eq.king.or..not.out1file)
673 & write (iout,*) "Reference temperature for weights calculation:",
675 call reada(weightcard,"D0CM",d0cm,3.78d0)
676 call reada(weightcard,"AKCM",akcm,15.1d0)
677 call reada(weightcard,"AKTH",akth,11.0d0)
678 call reada(weightcard,"AKCT",akct,12.0d0)
679 call reada(weightcard,"V1SS",v1ss,-1.08d0)
680 call reada(weightcard,"V2SS",v2ss,7.61d0)
681 call reada(weightcard,"V3SS",v3ss,13.7d0)
682 call reada(weightcard,"EBR",ebr,-5.50D0)
683 call reada(weightcard,"DTRISS",dtriss,1.0D0)
684 call reada(weightcard,"ATRISS",atriss,0.3D0)
685 call reada(weightcard,"BTRISS",btriss,0.02D0)
686 call reada(weightcard,"CTRISS",ctriss,1.0D0)
687 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
689 dyn_ss_mask(i)=.false.
693 dyn_ssbond_ij(i,j)=1.0d300
696 call reada(weightcard,"HT",Ht,0.0D0)
698 ss_depth=ebr/wsc-0.25*eps(1,1)
699 Ht=Ht/wsc-0.25*eps(1,1)
700 akcm=akcm*wstrain/wsc
701 akth=akth*wstrain/wsc
702 akct=akct*wstrain/wsc
703 v1ss=v1ss*wstrain/wsc
704 v2ss=v2ss*wstrain/wsc
705 v3ss=v3ss*wstrain/wsc
707 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
710 if(me.eq.king.or..not.out1file) then
711 write (iout,*) "Parameters of the SS-bond potential:"
712 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
714 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
715 write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
716 write (iout,*)" HT",Ht
717 print *,'indpdb=',indpdb,' pdbref=',pdbref
719 if (indpdb.gt.0 .or. pdbref) then
720 read(inp,'(a)') pdbfile
721 if(me.eq.king.or..not.out1file)
722 & write (iout,'(2a)') 'PDB data will be read from file ',
723 & pdbfile(:ilen(pdbfile))
724 open(ipdbin,file=pdbfile,status='old',err=33)
726 33 write (iout,'(a)') 'Error opening PDB file.'
729 c write (iout,*) 'Begin reading pdb data'
732 c write (iout,*) 'Finished reading pdb data'
734 if(me.eq.king.or..not.out1file)
735 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
736 & ' nstart_sup=',nstart_sup
738 itype_pdb(i)=itype(i)
742 nct=nstart_sup+nsup-1
743 call contact(.false.,ncont_ref,icont_ref,co)
746 if(me.eq.king.or..not.out1file)
747 & write(iout,*)'Adding sidechains'
751 if (iti.ne.10 .and. itype(i).ne.ntyp1) then
754 do while (fail.and.nsi.le.maxsi)
755 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
758 if(fail) write(iout,*)'Adding sidechain failed for res ',
759 & i,' after ',nsi,' trials'
764 if (indpdb.eq.0) then
765 C Read sequence if not taken from the pdb file.
767 c print *,'nres=',nres
768 if (iscode.gt.0) then
769 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
771 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
773 C Convert sequence to numeric code
775 itype(i)=rescode(i,sequence(i),iscode)
777 C Assign initial virtual bond lengths
783 vbld(i+nres)=dsc(iabs(itype(i)))
784 vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
785 c write (iout,*) "i",i," itype",itype(i),
786 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
790 c print '(20i4)',(itype(i),i=1,nres)
793 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
795 if (itype(i).eq.ntyp1) then
799 else if (iabs(itype(i+1)).ne.20) then
801 else if (iabs(itype(i)).ne.20) then
808 if(me.eq.king.or..not.out1file)then
809 write (iout,*) "ITEL"
811 write (iout,*) i,itype(i),itel(i)
813 print *,'Call Read_Bridge.'
816 C 8/13/98 Set limits to generating the dihedral angles
821 read (inp,*) ndih_constr
822 if (ndih_constr.gt.0) then
824 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
825 if(me.eq.king.or..not.out1file)then
827 & 'There are',ndih_constr,' constraints on phi angles.'
829 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
833 phi0(i)=deg2rad*phi0(i)
834 drange(i)=deg2rad*drange(i)
836 if(me.eq.king.or..not.out1file)
837 & write (iout,*) 'FTORS',ftors
840 phibound(1,ii) = phi0(i)-drange(i)
841 phibound(2,ii) = phi0(i)+drange(i)
848 write (iout,'(a)') 'Boundaries in phi angle sampling:'
850 write (iout,'(a3,i5,2f10.1)')
851 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
857 cd print *,'NNT=',NNT,' NCT=',NCT
858 if (itype(1).eq.ntyp1) nnt=2
859 if (itype(nres).eq.ntyp1) nct=nct-1
861 if(me.eq.king.or..not.out1file)
862 & write (iout,'(a,i3)') 'nsup=',nsup
864 if (nsup.le.(nct-nnt+1)) then
865 do i=0,nct-nnt+1-nsup
866 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
872 & 'Error - sequences to be superposed do not match.'
875 do i=0,nsup-(nct-nnt+1)
876 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
878 nstart_sup=nstart_sup+i
884 & 'Error - sequences to be superposed do not match.'
887 if (nsup.eq.0) nsup=nct-nnt
888 if (nstart_sup.eq.0) nstart_sup=nnt
889 if (nstart_seq.eq.0) nstart_seq=nnt
890 if(me.eq.king.or..not.out1file)
891 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
892 & ' nstart_seq=',nstart_seq
894 c--- Zscore rms -------
895 if (nz_start.eq.0) nz_start=nnt
896 if (nz_end.eq.0 .and. nsup.gt.0) then
898 else if (nz_end.eq.0) then
901 if(me.eq.king.or..not.out1file)then
902 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
903 write (iout,*) 'IZ_SC=',iz_sc
905 c----------------------
908 if (.not.pdbref) then
909 call read_angles(inp,*38)
911 38 write (iout,'(a)') 'Error reading reference structure.'
913 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
914 stop 'Error reading reference structure'
918 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
928 call contact(.true.,ncont_ref,icont_ref,co)
930 c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
932 if (constr_dist.gt.0) call read_dist_constr
933 write (iout,*) "After read_dist_constr nhpb",nhpb
935 if(me.eq.king.or..not.out1file)
936 & write (iout,*) 'Contact order:',co
938 if(me.eq.king.or..not.out1file)
939 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
942 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
944 if(me.eq.king.or..not.out1file)
945 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
946 & icont_ref(1,i),' ',
947 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
951 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
952 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
953 & modecalc.ne.10) then
954 C If input structure hasn't been supplied from the PDB file read or generate
956 if (iranconf.eq.0 .and. .not. extconf) then
957 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
958 & write (iout,'(a)') 'Initial geometry will be read in.'
960 read(inp,'(8f10.5)',end=36,err=36)
961 & ((c(l,k),l=1,3),k=1,nres),
962 & ((c(l,k+nres),l=1,3),k=nnt,nct)
963 write (iout,*) "Exit READ_CART"
964 write (iout,'(8f10.5)')
965 & ((c(l,k),l=1,3),k=1,nres),
966 & ((c(l,k+nres),l=1,3),k=nnt,nct)
967 call int_from_cart1(.true.)
968 write (iout,*) "Finish INT_TO_CART"
971 dc(j,i)=c(j,i+1)-c(j,i)
972 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
976 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
978 dc(j,i+nres)=c(j,i+nres)-c(j,i)
979 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
985 call read_angles(inp,*36)
988 36 write (iout,'(a)') 'Error reading angle file.'
990 call mpi_finalize( MPI_COMM_WORLD,IERR )
992 stop 'Error reading angle file.'
994 else if (extconf) then
995 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
996 & write (iout,'(a)') 'Extended chain initial geometry.'
998 theta(i)=90d0*deg2rad
1001 phi(i)=180d0*deg2rad
1004 alph(i)=110d0*deg2rad
1007 omeg(i)=-120d0*deg2rad
1008 if (itype(i).le.0) omeg(i)=-omeg(i)
1011 if(me.eq.king.or..not.out1file)
1012 & write (iout,'(a)') 'Random-generated initial geometry.'
1016 if (me.eq.king .or. fg_rank.eq.0 .and. (
1017 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1021 call gen_rand_conf(itmp,*30)
1023 30 write (iout,*) 'Failed to generate random conformation',
1024 & ', itrial=',itrial
1025 write (*,*) 'Processor:',me,
1026 & ' Failed to generate random conformation',
1036 write (iout,'(a,i3,a)') 'Processor:',me,
1037 & ' error in generating random conformation.'
1038 write (*,'(a,i3,a)') 'Processor:',me,
1039 & ' error in generating random conformation.'
1042 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1048 call gen_rand_conf(itmp,*30)
1050 30 write (iout,*) 'Failed to generate random conformation',
1051 & ', itrial=',itrial
1052 write (*,*) 'Failed to generate random conformation',
1053 & ', itrial=',itrial
1055 write (iout,'(a,i3,a)') 'Processor:',me,
1056 & ' error in generating random conformation.'
1057 write (*,'(a,i3,a)') 'Processor:',me,
1058 & ' error in generating random conformation.'
1063 elseif (modecalc.eq.4) then
1064 read (inp,'(a)') intinname
1065 open (intin,file=intinname,status='old',err=333)
1066 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1067 & write (iout,'(a)') 'intinname',intinname
1068 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1070 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1072 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1074 stop 'Error opening angle file.'
1078 C Generate distance constraints, if the PDB structure is to be regularized.
1079 if (nthread.gt.0) then
1080 call read_threadbase
1083 if (me.eq.king .or. .not. out1file)
1085 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1086 write (iout,'(/a,i3,a)')
1087 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1088 write (iout,'(20i4)') (iss(i),i=1,ns)
1090 write(iout,*)"Running with dynamic disulfide-bond formation"
1092 write (iout,'(/a/)') 'Pre-formed links are:'
1098 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1099 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1105 if (ns.gt.0.and.dyn_ss) then
1109 forcon(i-nss)=forcon(i)
1116 dyn_ss_mask(iss(i))=.true.
1119 if (i2ndstr.gt.0) call secstrp2dihc
1120 c call geom_to_var(nvar,x)
1121 c call etotal(energia(0))
1122 c call enerprint(energia(0))
1123 c call briefout(0,etot)
1125 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1126 cd write (iout,'(a)') 'Variable list:'
1127 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1129 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1130 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1131 & 'Processor',myrank,': end reading molecular data.'
1135 c--------------------------------------------------------------------------
1136 logical function seq_comp(itypea,itypeb,length)
1138 integer length,itypea(length),itypeb(length)
1141 if (itypea(i).ne.itypeb(i)) then
1149 c-----------------------------------------------------------------------------
1150 subroutine read_bridge
1151 C Read information about disulfide bridges.
1152 implicit real*8 (a-h,o-z)
1153 include 'DIMENSIONS'
1157 include 'COMMON.IOUNITS'
1158 include 'COMMON.GEO'
1159 include 'COMMON.VAR'
1160 include 'COMMON.INTERACT'
1161 include 'COMMON.LOCAL'
1162 include 'COMMON.NAMES'
1163 include 'COMMON.CHAIN'
1164 include 'COMMON.FFIELD'
1165 include 'COMMON.SBRIDGE'
1166 include 'COMMON.HEADER'
1167 include 'COMMON.CONTROL'
1168 include 'COMMON.DBASE'
1169 include 'COMMON.THREAD'
1170 include 'COMMON.TIME1'
1171 include 'COMMON.SETUP'
1172 C Read bridging residues.
1173 read (inp,*) ns,(iss(i),i=1,ns)
1175 if(me.eq.king.or..not.out1file)
1176 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1177 C Check whether the specified bridging residues are cystines.
1179 if (itype(iss(i)).ne.1) then
1180 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1181 & 'Do you REALLY think that the residue ',
1182 & restyp(itype(iss(i))),i,
1183 & ' can form a disulfide bridge?!!!'
1184 write (*,'(2a,i3,a)')
1185 & 'Do you REALLY think that the residue ',
1186 & restyp(itype(iss(i))),i,
1187 & ' can form a disulfide bridge?!!!'
1189 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1194 C Read preformed bridges.
1196 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1198 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1201 C Check if the residues involved in bridges are in the specified list of
1202 C bridging residues.
1205 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1206 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1207 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1208 & ' contains residues present in other pairs.'
1209 write (*,'(a,i3,a)') 'Disulfide pair',i,
1210 & ' contains residues present in other pairs.'
1212 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1218 if (ihpb(i).eq.iss(j)) goto 10
1220 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1223 if (jhpb(i).eq.iss(j)) goto 20
1225 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1231 ihpb(i)=ihpb(i)+nres
1232 jhpb(i)=jhpb(i)+nres
1238 c----------------------------------------------------------------------------
1239 subroutine read_x(kanal,*)
1240 implicit real*8 (a-h,o-z)
1241 include 'DIMENSIONS'
1242 include 'COMMON.GEO'
1243 include 'COMMON.VAR'
1244 include 'COMMON.CHAIN'
1245 include 'COMMON.IOUNITS'
1246 include 'COMMON.CONTROL'
1247 include 'COMMON.LOCAL'
1248 include 'COMMON.INTERACT'
1249 c Read coordinates from input
1251 read(kanal,'(8f10.5)',end=10,err=10)
1252 & ((c(l,k),l=1,3),k=1,nres),
1253 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1256 c(j,2*nres)=c(j,nres)
1258 call int_from_cart1(.false.)
1261 dc(j,i)=c(j,i+1)-c(j,i)
1262 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1266 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1268 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1269 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1277 c----------------------------------------------------------------------------
1278 subroutine read_threadbase
1279 implicit real*8 (a-h,o-z)
1280 include 'DIMENSIONS'
1281 include 'COMMON.IOUNITS'
1282 include 'COMMON.GEO'
1283 include 'COMMON.VAR'
1284 include 'COMMON.INTERACT'
1285 include 'COMMON.LOCAL'
1286 include 'COMMON.NAMES'
1287 include 'COMMON.CHAIN'
1288 include 'COMMON.FFIELD'
1289 include 'COMMON.SBRIDGE'
1290 include 'COMMON.HEADER'
1291 include 'COMMON.CONTROL'
1292 include 'COMMON.DBASE'
1293 include 'COMMON.THREAD'
1294 include 'COMMON.TIME1'
1295 C Read pattern database for threading.
1296 read (icbase,*) nseq
1298 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1299 & nres_base(2,i),nres_base(3,i)
1300 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1302 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1303 c & nres_base(2,i),nres_base(3,i)
1304 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1308 if (weidis.eq.0.0D0) weidis=0.1D0
1317 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1318 write (iout,'(a,i5)') 'nexcl: ',nexcl
1319 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1322 c------------------------------------------------------------------------------
1323 subroutine setup_var
1324 implicit real*8 (a-h,o-z)
1325 include 'DIMENSIONS'
1326 include 'COMMON.IOUNITS'
1327 include 'COMMON.GEO'
1328 include 'COMMON.VAR'
1329 include 'COMMON.INTERACT'
1330 include 'COMMON.LOCAL'
1331 include 'COMMON.NAMES'
1332 include 'COMMON.CHAIN'
1333 include 'COMMON.FFIELD'
1334 include 'COMMON.SBRIDGE'
1335 include 'COMMON.HEADER'
1336 include 'COMMON.CONTROL'
1337 include 'COMMON.DBASE'
1338 include 'COMMON.THREAD'
1339 include 'COMMON.TIME1'
1340 C Set up variable list.
1346 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1348 ialph(i,1)=nvar+nside
1352 if (indphi.gt.0) then
1354 else if (indback.gt.0) then
1359 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1362 c----------------------------------------------------------------------------
1363 subroutine gen_dist_constr
1364 C Generate CA distance constraints.
1365 implicit real*8 (a-h,o-z)
1366 include 'DIMENSIONS'
1367 include 'COMMON.IOUNITS'
1368 include 'COMMON.GEO'
1369 include 'COMMON.VAR'
1370 include 'COMMON.INTERACT'
1371 include 'COMMON.LOCAL'
1372 include 'COMMON.NAMES'
1373 include 'COMMON.CHAIN'
1374 include 'COMMON.FFIELD'
1375 include 'COMMON.SBRIDGE'
1376 include 'COMMON.HEADER'
1377 include 'COMMON.CONTROL'
1378 include 'COMMON.DBASE'
1379 include 'COMMON.THREAD'
1380 include 'COMMON.TIME1'
1381 dimension itype_pdb(maxres)
1382 common /pizda/ itype_pdb
1384 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1385 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1386 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1388 do i=nstart_sup,nstart_sup+nsup-1
1389 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1390 cd & ' seq_pdb', restyp(itype_pdb(i))
1391 do j=i+2,nstart_sup+nsup-1
1393 ihpb(nhpb)=i+nstart_seq-nstart_sup
1394 jhpb(nhpb)=j+nstart_seq-nstart_sup
1396 dhpb(nhpb)=dist(i,j)
1399 cd write (iout,'(a)') 'Distance constraints:'
1404 cd if (ii.gt.nres) then
1409 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1410 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1411 cd & dhpb(i),forcon(i)
1415 c----------------------------------------------------------------------------
1417 implicit real*8 (a-h,o-z)
1418 include 'DIMENSIONS'
1419 include 'COMMON.MAP'
1420 include 'COMMON.IOUNITS'
1421 character*3 angid(4) /'THE','PHI','ALP','OME'/
1422 character*80 mapcard,ucase
1424 read (inp,'(a)') mapcard
1425 mapcard=ucase(mapcard)
1426 if (index(mapcard,'PHI').gt.0) then
1428 else if (index(mapcard,'THE').gt.0) then
1430 else if (index(mapcard,'ALP').gt.0) then
1432 else if (index(mapcard,'OME').gt.0) then
1435 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1436 stop 'Error - illegal variable spec in MAP card.'
1438 call readi (mapcard,'RES1',res1(imap),0)
1439 call readi (mapcard,'RES2',res2(imap),0)
1440 if (res1(imap).eq.0) then
1441 res1(imap)=res2(imap)
1442 else if (res2(imap).eq.0) then
1443 res2(imap)=res1(imap)
1445 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1447 & 'Error - illegal definition of variable group in MAP.'
1448 stop 'Error - illegal definition of variable group in MAP.'
1450 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1451 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1452 call readi(mapcard,'NSTEP',nstep(imap),0)
1453 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1455 & 'Illegal boundary and/or step size specification in MAP.'
1456 stop 'Illegal boundary and/or step size specification in MAP.'
1461 c----------------------------------------------------------------------------
1463 implicit real*8 (a-h,o-z)
1464 include 'DIMENSIONS'
1465 include 'COMMON.IOUNITS'
1466 include 'COMMON.GEO'
1467 include 'COMMON.CSA'
1468 include 'COMMON.BANK'
1469 include 'COMMON.CONTROL'
1471 character*620 mcmcard
1472 call card_concat(mcmcard)
1474 call readi(mcmcard,'NCONF',nconf,50)
1475 call readi(mcmcard,'NADD',nadd,0)
1476 call readi(mcmcard,'JSTART',jstart,1)
1477 call readi(mcmcard,'JEND',jend,1)
1478 call readi(mcmcard,'NSTMAX',nstmax,500000)
1479 call readi(mcmcard,'N0',n0,1)
1480 call readi(mcmcard,'N1',n1,6)
1481 call readi(mcmcard,'N2',n2,4)
1482 call readi(mcmcard,'N3',n3,0)
1483 call readi(mcmcard,'N4',n4,0)
1484 call readi(mcmcard,'N5',n5,0)
1485 call readi(mcmcard,'N6',n6,10)
1486 call readi(mcmcard,'N7',n7,0)
1487 call readi(mcmcard,'N8',n8,0)
1488 call readi(mcmcard,'N9',n9,0)
1489 call readi(mcmcard,'N14',n14,0)
1490 call readi(mcmcard,'N15',n15,0)
1491 call readi(mcmcard,'N16',n16,0)
1492 call readi(mcmcard,'N17',n17,0)
1493 call readi(mcmcard,'N18',n18,0)
1495 vdisulf=(index(mcmcard,'DYNSS').gt.0)
1497 call readi(mcmcard,'NDIFF',ndiff,2)
1498 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1499 call readi(mcmcard,'IS1',is1,1)
1500 call readi(mcmcard,'IS2',is2,8)
1501 call readi(mcmcard,'NRAN0',nran0,4)
1502 call readi(mcmcard,'NRAN1',nran1,2)
1503 call readi(mcmcard,'IRR',irr,1)
1504 call readi(mcmcard,'NSEED',nseed,20)
1505 call readi(mcmcard,'NTOTAL',ntotal,10000)
1506 call reada(mcmcard,'CUT1',cut1,2.0d0)
1507 call reada(mcmcard,'CUT2',cut2,5.0d0)
1508 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1509 call readi(mcmcard,'ICMAX',icmax,3)
1510 call readi(mcmcard,'IRESTART',irestart,0)
1511 c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1514 call reada(mcmcard,'DELE',dele,20.0d0)
1515 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1516 call readi(mcmcard,'IREF',iref,0)
1517 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1518 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1519 call readi(mcmcard,'NCONF_IN',nconf_in,0)
1520 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1521 write (iout,*) "NCONF_IN",nconf_in
1524 c----------------------------------------------------------------------------
1525 cfmc subroutine mcmfread
1526 cfmc implicit real*8 (a-h,o-z)
1527 cfmc include 'DIMENSIONS'
1528 cfmc include 'COMMON.MCMF'
1529 cfmc include 'COMMON.IOUNITS'
1530 cfmc include 'COMMON.GEO'
1531 cfmc character*80 ucase
1532 cfmc character*620 mcmcard
1533 cfmc call card_concat(mcmcard)
1535 cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
1536 cfmc write(iout,*)'MAXRANT=',maxrant
1537 cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1538 cfmc write(iout,*)'MAXFAM=',maxfam
1539 cfmc call readi(mcmcard,'NNET1',nnet1,5)
1540 cfmc write(iout,*)'NNET1=',nnet1
1541 cfmc call readi(mcmcard,'NNET2',nnet2,4)
1542 cfmc write(iout,*)'NNET2=',nnet2
1543 cfmc call readi(mcmcard,'NNET3',nnet3,4)
1544 cfmc write(iout,*)'NNET3=',nnet3
1545 cfmc call readi(mcmcard,'ILASTT',ilastt,0)
1546 cfmc write(iout,*)'ILASTT=',ilastt
1547 cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1548 cfmc write(iout,*)'MAXSTR=',maxstr
1549 cfmc maxstr_f=maxstr/maxfam
1550 cfmc write(iout,*)'MAXSTR_F=',maxstr_f
1551 cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
1552 cfmc write(iout,*)'NMCMF=',nmcmf
1553 cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1554 cfmc write(iout,*)'IFOCUS=',ifocus
1555 cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1556 cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
1557 cfmc call readi(mcmcard,'INTPRT',intprt,1000)
1558 cfmc write(iout,*)'INTPRT=',intprt
1559 cfmc call readi(mcmcard,'IPRT',iprt,100)
1560 cfmc write(iout,*)'IPRT=',iprt
1561 cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
1562 cfmc write(iout,*)'IMAXTR=',imaxtr
1563 cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
1564 cfmc write(iout,*)'MAXEVEN=',maxeven
1565 cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1566 cfmc write(iout,*)'MAXEVEN1=',maxeven1
1567 cfmc call readi(mcmcard,'INIMIN',inimin,200)
1568 cfmc write(iout,*)'INIMIN=',inimin
1569 cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1570 cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
1571 cfmc call readi(mcmcard,'NTHREAD',nthread,5)
1572 cfmc write(iout,*)'NTHREAD=',nthread
1573 cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1574 cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1575 cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
1576 cfmc write(iout,*)'MAXPERT=',maxpert
1577 cfmc call readi(mcmcard,'IRMSD',irmsd,1)
1578 cfmc write(iout,*)'IRMSD=',irmsd
1579 cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1580 cfmc write(iout,*)'DENEMIN=',denemin
1581 cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1582 cfmc write(iout,*)'RCUT1S=',rcut1s
1583 cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1584 cfmc write(iout,*)'RCUT1E=',rcut1e
1585 cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1586 cfmc write(iout,*)'RCUT2S=',rcut2s
1587 cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1588 cfmc write(iout,*)'RCUT2E=',rcut2e
1589 cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1590 cfmc write(iout,*)'DPERT1=',d_pert1
1591 cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1592 cfmc write(iout,*)'DPERT1A=',d_pert1a
1593 cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1594 cfmc write(iout,*)'DPERT2=',d_pert2
1595 cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1596 cfmc write(iout,*)'DPERT2A=',d_pert2a
1597 cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1598 cfmc write(iout,*)'DPERT2B=',d_pert2b
1599 cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1600 cfmc write(iout,*)'DPERT2C=',d_pert2c
1601 cfmc d_pert1=deg2rad*d_pert1
1602 cfmc d_pert1a=deg2rad*d_pert1a
1603 cfmc d_pert2=deg2rad*d_pert2
1604 cfmc d_pert2a=deg2rad*d_pert2a
1605 cfmc d_pert2b=deg2rad*d_pert2b
1606 cfmc d_pert2c=deg2rad*d_pert2c
1607 cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1608 cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
1609 cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1610 cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
1611 cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1612 cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1613 cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1614 cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1615 cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1616 cfmc write(iout,*)'RCUTINI=',rcutini
1617 cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
1618 cfmc write(iout,*)'GRAT=',grat
1619 cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1620 cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
1624 c----------------------------------------------------------------------------
1626 implicit real*8 (a-h,o-z)
1627 include 'DIMENSIONS'
1628 include 'COMMON.MCM'
1629 include 'COMMON.MCE'
1630 include 'COMMON.IOUNITS'
1632 character*320 mcmcard
1633 call card_concat(mcmcard)
1634 call readi(mcmcard,'MAXACC',maxacc,100)
1635 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1636 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1637 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1638 call readi(mcmcard,'MAXREPM',maxrepm,200)
1639 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1640 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1641 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1642 call reada(mcmcard,'E_UP',e_up,5.0D0)
1643 call reada(mcmcard,'DELTE',delte,0.1D0)
1644 call readi(mcmcard,'NSWEEP',nsweep,5)
1645 call readi(mcmcard,'NSTEPH',nsteph,0)
1646 call readi(mcmcard,'NSTEPC',nstepc,0)
1647 call reada(mcmcard,'TMIN',tmin,298.0D0)
1648 call reada(mcmcard,'TMAX',tmax,298.0D0)
1649 call readi(mcmcard,'NWINDOW',nwindow,0)
1650 call readi(mcmcard,'PRINT_MC',print_mc,0)
1651 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1652 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1653 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1654 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1655 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1656 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1657 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1658 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1659 if (nwindow.gt.0) then
1660 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1662 winlen(i)=winend(i)-winstart(i)+1
1665 if (tmax.lt.tmin) tmax=tmin
1666 if (tmax.eq.tmin) then
1670 if (nstepc.gt.0 .and. nsteph.gt.0) then
1671 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1672 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1674 C Probabilities of different move types
1675 sumpro_type(0)=0.0D0
1676 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1677 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1678 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1679 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1680 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1681 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1682 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1684 print *,'i',i,' sumprotype',sumpro_type(i)
1685 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1686 print *,'i',i,' sumprotype',sumpro_type(i)
1690 c----------------------------------------------------------------------------
1691 subroutine read_minim
1692 implicit real*8 (a-h,o-z)
1693 include 'DIMENSIONS'
1694 include 'COMMON.MINIM'
1695 include 'COMMON.IOUNITS'
1697 character*320 minimcard
1698 call card_concat(minimcard)
1699 call readi(minimcard,'MAXMIN',maxmin,2000)
1700 call readi(minimcard,'MAXFUN',maxfun,5000)
1701 call readi(minimcard,'MINMIN',minmin,maxmin)
1702 call readi(minimcard,'MINFUN',minfun,maxmin)
1703 call reada(minimcard,'TOLF',tolf,1.0D-2)
1704 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1705 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1706 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1707 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1708 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1709 & 'Options in energy minimization:'
1710 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1711 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1712 & 'MinMin:',MinMin,' MinFun:',MinFun,
1713 & ' TolF:',TolF,' RTolF:',RTolF
1716 c----------------------------------------------------------------------------
1717 subroutine read_angles(kanal,*)
1718 implicit real*8 (a-h,o-z)
1719 include 'DIMENSIONS'
1720 include 'COMMON.GEO'
1721 include 'COMMON.VAR'
1722 include 'COMMON.CHAIN'
1723 include 'COMMON.IOUNITS'
1724 include 'COMMON.CONTROL'
1725 c Read angles from input
1727 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1728 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1729 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1730 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1733 c 9/7/01 avoid 180 deg valence angle
1734 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1736 theta(i)=deg2rad*theta(i)
1737 phi(i)=deg2rad*phi(i)
1738 alph(i)=deg2rad*alph(i)
1739 omeg(i)=deg2rad*omeg(i)
1744 c----------------------------------------------------------------------------
1745 subroutine reada(rekord,lancuch,wartosc,default)
1747 character*(*) rekord,lancuch
1748 double precision wartosc,default
1751 iread=index(rekord,lancuch)
1752 if (iread.eq.0) then
1756 iread=iread+ilen(lancuch)+1
1757 read (rekord(iread:),*,err=10,end=10) wartosc
1762 c----------------------------------------------------------------------------
1763 subroutine readi(rekord,lancuch,wartosc,default)
1765 character*(*) rekord,lancuch
1766 integer wartosc,default
1769 iread=index(rekord,lancuch)
1770 if (iread.eq.0) then
1774 iread=iread+ilen(lancuch)+1
1775 read (rekord(iread:),*,err=10,end=10) wartosc
1780 c----------------------------------------------------------------------------
1781 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1784 integer tablica(dim),default
1785 character*(*) rekord,lancuch
1792 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1793 if (iread.eq.0) return
1794 iread=iread+ilen(lancuch)+1
1795 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1798 c----------------------------------------------------------------------------
1799 subroutine multreada(rekord,lancuch,tablica,dim,default)
1802 double precision tablica(dim),default
1803 character*(*) rekord,lancuch
1810 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1811 if (iread.eq.0) return
1812 iread=iread+ilen(lancuch)+1
1813 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1816 c----------------------------------------------------------------------------
1817 subroutine openunits
1818 implicit real*8 (a-h,o-z)
1819 include 'DIMENSIONS'
1822 character*16 form,nodename
1825 include 'COMMON.SETUP'
1826 include 'COMMON.IOUNITS'
1828 include 'COMMON.CONTROL'
1829 integer lenpre,lenpot,ilen,lentmp
1831 character*3 out1file_text,ucase
1834 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1835 call getenv_loc("PREFIX",prefix)
1837 call getenv_loc("POT",pot)
1838 call getenv_loc("DIRTMP",tmpdir)
1839 call getenv_loc("CURDIR",curdir)
1840 call getenv_loc("OUT1FILE",out1file_text)
1841 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1842 out1file_text=ucase(out1file_text)
1843 if (out1file_text(1:1).eq."Y") then
1846 out1file=fg_rank.gt.0
1851 if (lentmp.gt.0) then
1852 write (*,'(80(1h!))')
1853 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
1854 write (*,'(80(1h!))')
1855 write (*,*)"All output files will be on node /tmp directory."
1857 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1858 if (me.eq.king) then
1859 write (*,*) "The master node is ",nodename
1860 else if (fg_rank.eq.0) then
1861 write (*,*) "I am the CG slave node ",nodename
1863 write (*,*) "I am the FG slave node ",nodename
1866 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1867 lenpre = lentmp+lenpre+1
1869 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
1870 C Get the names and open the input files
1871 #if defined(WINIFL) || defined(WINPGI)
1872 open(1,file=pref_orig(:ilen(pref_orig))//
1873 & '.inp',status='old',readonly,shared)
1874 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1875 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1876 C Get parameter filenames and open the parameter files.
1877 call getenv_loc('BONDPAR',bondname)
1878 open (ibond,file=bondname,status='old',readonly,shared)
1879 call getenv_loc('THETPAR',thetname)
1880 open (ithep,file=thetname,status='old',readonly,shared)
1881 call getenv_loc('ROTPAR',rotname)
1882 open (irotam,file=rotname,status='old',readonly,shared)
1883 call getenv_loc('TORPAR',torname)
1884 open (itorp,file=torname,status='old',readonly,shared)
1885 call getenv_loc('TORDPAR',tordname)
1886 open (itordp,file=tordname,status='old',readonly,shared)
1887 call getenv_loc('FOURIER',fouriername)
1888 open (ifourier,file=fouriername,status='old',readonly,shared)
1889 call getenv_loc('ELEPAR',elename)
1890 open (ielep,file=elename,status='old',readonly,shared)
1891 call getenv_loc('SIDEPAR',sidename)
1892 open (isidep,file=sidename,status='old',readonly,shared)
1893 #elif (defined CRAY) || (defined AIX)
1894 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1896 c print *,"Processor",myrank," opened file 1"
1897 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1898 c print *,"Processor",myrank," opened file 9"
1899 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1900 C Get parameter filenames and open the parameter files.
1901 call getenv_loc('BONDPAR',bondname)
1902 open (ibond,file=bondname,status='old',action='read')
1903 c print *,"Processor",myrank," opened file IBOND"
1904 call getenv_loc('THETPAR',thetname)
1905 open (ithep,file=thetname,status='old',action='read')
1906 c print *,"Processor",myrank," opened file ITHEP"
1907 call getenv_loc('ROTPAR',rotname)
1908 open (irotam,file=rotname,status='old',action='read')
1909 c print *,"Processor",myrank," opened file IROTAM"
1910 call getenv_loc('TORPAR',torname)
1911 open (itorp,file=torname,status='old',action='read')
1912 c print *,"Processor",myrank," opened file ITORP"
1913 call getenv_loc('TORDPAR',tordname)
1914 open (itordp,file=tordname,status='old',action='read')
1915 c print *,"Processor",myrank," opened file ITORDP"
1916 call getenv_loc('SCCORPAR',sccorname)
1917 open (isccor,file=sccorname,status='old',action='read')
1918 c print *,"Processor",myrank," opened file ISCCOR"
1919 call getenv_loc('FOURIER',fouriername)
1920 open (ifourier,file=fouriername,status='old',action='read')
1921 c print *,"Processor",myrank," opened file IFOURIER"
1922 call getenv_loc('ELEPAR',elename)
1923 open (ielep,file=elename,status='old',action='read')
1924 c print *,"Processor",myrank," opened file IELEP"
1925 call getenv_loc('SIDEPAR',sidename)
1926 open (isidep,file=sidename,status='old',action='read')
1927 c print *,"Processor",myrank," opened file ISIDEP"
1928 c print *,"Processor",myrank," opened parameter files"
1930 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
1931 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1932 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1933 C Get parameter filenames and open the parameter files.
1934 call getenv_loc('BONDPAR',bondname)
1935 open (ibond,file=bondname,status='old')
1936 call getenv_loc('THETPAR',thetname)
1937 open (ithep,file=thetname,status='old')
1938 call getenv_loc('ROTPAR',rotname)
1939 open (irotam,file=rotname,status='old')
1940 call getenv_loc('TORPAR',torname)
1941 open (itorp,file=torname,status='old')
1942 call getenv_loc('TORDPAR',tordname)
1943 open (itordp,file=tordname,status='old')
1944 call getenv_loc('SCCORPAR',sccorname)
1945 open (isccor,file=sccorname,status='old')
1946 call getenv_loc('FOURIER',fouriername)
1947 open (ifourier,file=fouriername,status='old')
1948 call getenv_loc('ELEPAR',elename)
1949 open (ielep,file=elename,status='old')
1950 call getenv_loc('SIDEPAR',sidename)
1951 open (isidep,file=sidename,status='old')
1953 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1955 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1956 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1957 C Get parameter filenames and open the parameter files.
1958 call getenv_loc('BONDPAR',bondname)
1959 open (ibond,file=bondname,status='old',readonly)
1960 call getenv_loc('THETPAR',thetname)
1961 open (ithep,file=thetname,status='old',readonly)
1962 call getenv_loc('ROTPAR',rotname)
1963 open (irotam,file=rotname,status='old',readonly)
1964 call getenv_loc('TORPAR',torname)
1965 open (itorp,file=torname,status='old',readonly)
1966 call getenv_loc('TORDPAR',tordname)
1967 open (itordp,file=tordname,status='old',readonly)
1968 call getenv_loc('SCCORPAR',sccorname)
1969 open (isccor,file=sccorname,status='old',readonly)
1971 call getenv_loc('THETPARPDB',thetname_pdb)
1972 print *,"thetname_pdb ",thetname_pdb
1973 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
1974 print *,ithep_pdb," opened"
1976 call getenv_loc('FOURIER',fouriername)
1977 open (ifourier,file=fouriername,status='old',readonly)
1978 call getenv_loc('ELEPAR',elename)
1979 open (ielep,file=elename,status='old',readonly)
1980 call getenv_loc('SIDEPAR',sidename)
1981 open (isidep,file=sidename,status='old',readonly)
1983 call getenv_loc('ROTPARPDB',rotname_pdb)
1984 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
1989 C 8/9/01 In the newest version SCp interaction constants are read from a file
1990 C Use -DOLDSCP to use hard-coded constants instead.
1992 call getenv_loc('SCPPAR',scpname)
1993 #if defined(WINIFL) || defined(WINPGI)
1994 open (iscpp,file=scpname,status='old',readonly,shared)
1995 #elif (defined CRAY) || (defined AIX)
1996 open (iscpp,file=scpname,status='old',action='read')
1998 open (iscpp,file=scpname,status='old')
2000 open (iscpp,file=scpname,status='old',readonly)
2003 call getenv_loc('PATTERN',patname)
2004 #if defined(WINIFL) || defined(WINPGI)
2005 open (icbase,file=patname,status='old',readonly,shared)
2006 #elif (defined CRAY) || (defined AIX)
2007 open (icbase,file=patname,status='old',action='read')
2009 open (icbase,file=patname,status='old')
2011 open (icbase,file=patname,status='old',readonly)
2014 C Open output file only for CG processes
2015 c print *,"Processor",myrank," fg_rank",fg_rank
2016 if (fg_rank.eq.0) then
2018 if (nodes.eq.1) then
2021 npos = dlog10(dfloat(nodes-1))+1
2023 if (npos.lt.3) npos=3
2024 write (liczba,'(i1)') npos
2025 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2027 write (liczba,form) me
2028 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2029 & liczba(:ilen(liczba))
2030 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2032 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2034 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2035 & liczba(:ilen(liczba))//'.mol2'
2036 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2037 & liczba(:ilen(liczba))//'.stat'
2039 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2040 & //liczba(:ilen(liczba))//'.stat')
2041 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2044 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2045 & liczba(:ilen(liczba))//'.const'
2050 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2051 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2052 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2053 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2054 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2056 & call copy_to_tmp(pref_orig(:ile(pref_orig))//'_'//pot(:lenpot)//
2058 rest2name=prefix(:ilen(prefix))//'.rst'
2060 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2063 #if defined(AIX) || defined(PGI)
2064 if (me.eq.king .or. .not. out1file)
2065 & open(iout,file=outname,status='unknown')
2067 if (fg_rank.gt.0) then
2068 write (liczba,'(i3.3)') myrank/nfgtasks
2069 write (ll,'(bz,i3.3)') fg_rank
2070 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2075 open(igeom,file=intname,status='unknown',position='append')
2076 open(ipdb,file=pdbname,status='unknown')
2077 open(imol2,file=mol2name,status='unknown')
2078 open(istat,file=statname,status='unknown',position='append')
2080 c1out open(iout,file=outname,status='unknown')
2083 if (me.eq.king .or. .not.out1file)
2084 & open(iout,file=outname,status='unknown')
2086 if (fg_rank.gt.0) then
2087 write (liczba,'(i3.3)') myrank/nfgtasks
2088 write (ll,'(bz,i3.3)') fg_rank
2089 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2094 open(igeom,file=intname,status='unknown',access='append')
2095 open(ipdb,file=pdbname,status='unknown')
2096 open(imol2,file=mol2name,status='unknown')
2097 open(istat,file=statname,status='unknown',access='append')
2099 c1out open(iout,file=outname,status='unknown')
2102 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2103 csa_seed=prefix(:lenpre)//'.CSA.seed'
2104 csa_history=prefix(:lenpre)//'.CSA.history'
2105 csa_bank=prefix(:lenpre)//'.CSA.bank'
2106 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2107 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2108 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2109 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2110 csa_int=prefix(:lenpre)//'.int'
2111 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2112 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2113 csa_in=prefix(:lenpre)//'.CSA.in'
2114 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2117 write (iout,'(80(1h-))')
2118 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2119 write (iout,'(80(1h-))')
2120 write (iout,*) "Input file : ",
2121 & pref_orig(:ilen(pref_orig))//'.inp'
2122 write (iout,*) "Output file : ",
2123 & outname(:ilen(outname))
2125 write (iout,*) "Sidechain potential file : ",
2126 & sidename(:ilen(sidename))
2128 write (iout,*) "SCp potential file : ",
2129 & scpname(:ilen(scpname))
2131 write (iout,*) "Electrostatic potential file : ",
2132 & elename(:ilen(elename))
2133 write (iout,*) "Cumulant coefficient file : ",
2134 & fouriername(:ilen(fouriername))
2135 write (iout,*) "Torsional parameter file : ",
2136 & torname(:ilen(torname))
2137 write (iout,*) "Double torsional parameter file : ",
2138 & tordname(:ilen(tordname))
2139 write (iout,*) "SCCOR parameter file : ",
2140 & sccorname(:ilen(sccorname))
2141 write (iout,*) "Bond & inertia constant file : ",
2142 & bondname(:ilen(bondname))
2143 write (iout,*) "Bending parameter file : ",
2144 & thetname(:ilen(thetname))
2145 write (iout,*) "Rotamer parameter file : ",
2146 & rotname(:ilen(rotname))
2147 write (iout,*) "Thetpdb parameter file : ",
2148 & thetname_pdb(:ilen(thetname_pdb))
2149 write (iout,*) "Threading database : ",
2150 & patname(:ilen(patname))
2152 &write (iout,*)" DIRTMP : ",
2154 write (iout,'(80(1h-))')
2158 c----------------------------------------------------------------------------
2159 subroutine card_concat(card)
2160 implicit real*8 (a-h,o-z)
2161 include 'DIMENSIONS'
2162 include 'COMMON.IOUNITS'
2164 character*80 karta,ucase
2166 read (inp,'(a)') karta
2169 do while (karta(80:80).eq.'&')
2170 card=card(:ilen(card)+1)//karta(:79)
2171 read (inp,'(a)') karta
2174 card=card(:ilen(card)+1)//karta
2177 c----------------------------------------------------------------------------------
2179 implicit real*8 (a-h,o-z)
2180 include 'DIMENSIONS'
2181 include 'COMMON.CHAIN'
2182 include 'COMMON.IOUNITS'
2184 open(irest2,file=rest2name,status='unknown')
2185 read(irest2,*) totT,EK,potE,totE,t_bath
2187 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2190 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2193 read (irest2,*) iset
2198 c---------------------------------------------------------------------------------
2199 subroutine read_fragments
2200 implicit real*8 (a-h,o-z)
2201 include 'DIMENSIONS'
2205 include 'COMMON.SETUP'
2206 include 'COMMON.CHAIN'
2207 include 'COMMON.IOUNITS'
2209 include 'COMMON.CONTROL'
2210 read(inp,*) nset,nfrag,npair,nfrag_back
2211 if(me.eq.king.or..not.out1file)
2212 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2213 & " nfrag_back",nfrag_back
2215 read(inp,*) mset(iset)
2217 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2219 if(me.eq.king.or..not.out1file)
2220 & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2221 & ifrag(2,i,iset), qinfrag(i,iset)
2224 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2226 if(me.eq.king.or..not.out1file)
2227 & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2228 & ipair(2,i,iset), qinpair(i,iset)
2231 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2232 & wfrag_back(3,i,iset),
2233 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2234 if(me.eq.king.or..not.out1file)
2235 & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2236 & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2241 c-------------------------------------------------------------------------------
2242 subroutine read_dist_constr
2243 implicit real*8 (a-h,o-z)
2244 include 'DIMENSIONS'
2248 include 'COMMON.SETUP'
2249 include 'COMMON.CONTROL'
2250 include 'COMMON.CHAIN'
2251 include 'COMMON.IOUNITS'
2252 include 'COMMON.SBRIDGE'
2253 integer ifrag_(2,100),ipair_(2,100)
2254 double precision wfrag_(100),wpair_(100)
2255 character*500 controlcard
2256 c write (iout,*) "Calling read_dist_constr"
2257 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2259 call card_concat(controlcard)
2260 call readi(controlcard,"NFRAG",nfrag_,0)
2261 call readi(controlcard,"NPAIR",npair_,0)
2262 call readi(controlcard,"NDIST",ndist_,0)
2263 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2264 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2265 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2266 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2267 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2268 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2269 c write (iout,*) "IFRAG"
2271 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2273 c write (iout,*) "IPAIR"
2275 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2279 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2280 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2281 & ifrag_(2,i)=nstart_sup+nsup-1
2282 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2284 if (wfrag_(i).gt.0.0d0) then
2285 do j=ifrag_(1,i),ifrag_(2,i)-1
2286 do k=j+1,ifrag_(2,i)
2287 c write (iout,*) "j",j," k",k
2289 if (constr_dist.eq.1) then
2294 forcon(nhpb)=wfrag_(i)
2295 else if (constr_dist.eq.2) then
2296 if (ddjk.le.dist_cut) then
2301 forcon(nhpb)=wfrag_(i)
2308 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2311 if (.not.out1file .or. me.eq.king)
2312 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2313 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2315 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2316 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2323 if (wpair_(i).gt.0.0d0) then
2331 do j=ifrag_(1,ii),ifrag_(2,ii)
2332 do k=ifrag_(1,jj),ifrag_(2,jj)
2336 forcon(nhpb)=wpair_(i)
2337 dhpb(nhpb)=dist(j,k)
2339 if (.not.out1file .or. me.eq.king)
2340 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2341 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2343 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2344 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2351 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2352 if (forcon(nhpb+1).gt.0.0d0) then
2354 dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2356 if (.not.out1file .or. me.eq.king)
2357 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2358 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2360 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2361 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2368 c-------------------------------------------------------------------------------
2370 subroutine flush(iu)
2375 subroutine flush(iu)
2380 c------------------------------------------------------------------------------
2381 subroutine copy_to_tmp(source)
2382 include "DIMENSIONS"
2383 include "COMMON.IOUNITS"
2384 character*(*) source
2385 character* 256 tmpfile
2389 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
2390 inquire(file=tmpfile,exist=ex)
2392 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
2393 & " to temporary directory..."
2394 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
2395 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
2399 c------------------------------------------------------------------------------
2400 subroutine move_from_tmp(source)
2401 include "DIMENSIONS"
2402 include "COMMON.IOUNITS"
2403 character*(*) source
2406 write (*,*) "Moving ",source(:ilen(source)),
2407 & " from temporary directory to working directory"
2408 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
2409 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
2412 c------------------------------------------------------------------------------
2413 subroutine random_init(seed)
2415 C Initialize random number generator
2417 implicit real*8 (a-h,o-z)
2418 include 'DIMENSIONS'
2421 logical OKRandom, prng_restart
2423 integer iseed_array(4)
2425 include 'COMMON.IOUNITS'
2426 include 'COMMON.TIME1'
2427 include 'COMMON.THREAD'
2428 include 'COMMON.SBRIDGE'
2429 include 'COMMON.CONTROL'
2430 include 'COMMON.MCM'
2431 include 'COMMON.MAP'
2432 include 'COMMON.HEADER'
2433 include 'COMMON.CSA'
2434 include 'COMMON.CHAIN'
2435 include 'COMMON.MUCA'
2437 include 'COMMON.FFIELD'
2438 include 'COMMON.SETUP'
2439 iseed=-dint(dabs(seed))
2440 if (iseed.eq.0) then
2441 write (iout,'(/80(1h*)/20x,a/80(1h*))')
2442 & 'Random seed undefined. The program will stop.'
2443 write (*,'(/80(1h*)/20x,a/80(1h*))')
2444 & 'Random seed undefined. The program will stop.'
2446 call mpi_finalize(mpi_comm_world,ierr)
2448 stop 'Bad random seed.'
2451 if (fg_rank.eq.0) then
2455 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
2456 OKRandom = prng_restart(me,iseed)
2459 tmp=65536.0d0**(4-i)
2460 iseed_array(i) = dint(seed/tmp)
2461 seed=seed-iseed_array(i)*tmp
2464 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
2465 & (iseed_array(i),i=1,4)
2466 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
2467 & (iseed_array(i),i=1,4)
2468 OKRandom = prng_restart(me,iseed_array)
2471 c r1 = prng_next(me)
2472 r1=ran_number(0.0D0,1.0D0)
2474 & write (iout,*) 'ran_num',r1
2475 if (r1.lt.0.0d0) OKRandom=.false.
2477 if (.not.OKRandom) then
2478 write (iout,*) 'PRNG IS NOT WORKING!!!'
2479 print *,'PRNG IS NOT WORKING!!!'
2482 call mpi_abort(mpi_comm_world,error_msg,ierr)
2485 write (iout,*) 'too many processors for parallel prng'
2486 write (*,*) 'too many processors for parallel prng'
2494 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)