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 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
685 dyn_ss_mask(i)=.false.
689 dyn_ssbond_ij(i,j)=1.0d300
692 call reada(weightcard,"HT",Ht,0.0D0)
694 ss_depth=ebr/wsc-0.25*eps(1,1)
695 Ht=Ht/wsc-0.25*eps(1,1)
696 akcm=akcm*wstrain/wsc
697 akth=akth*wstrain/wsc
698 akct=akct*wstrain/wsc
699 v1ss=v1ss*wstrain/wsc
700 v2ss=v2ss*wstrain/wsc
701 v3ss=v3ss*wstrain/wsc
703 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
706 if(me.eq.king.or..not.out1file) then
707 write (iout,*) "Parameters of the SS-bond potential:"
708 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
710 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
711 write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
712 write (iout,*)" HT",Ht
713 print *,'indpdb=',indpdb,' pdbref=',pdbref
715 if (indpdb.gt.0 .or. pdbref) then
716 read(inp,'(a)') pdbfile
717 if(me.eq.king.or..not.out1file)
718 & write (iout,'(2a)') 'PDB data will be read from file ',
719 & pdbfile(:ilen(pdbfile))
720 open(ipdbin,file=pdbfile,status='old',err=33)
722 33 write (iout,'(a)') 'Error opening PDB file.'
725 c print *,'Begin reading pdb data'
727 c print *,'Finished reading pdb data'
728 if(me.eq.king.or..not.out1file)
729 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
730 & ' nstart_sup=',nstart_sup
732 itype_pdb(i)=itype(i)
736 nct=nstart_sup+nsup-1
737 call contact(.false.,ncont_ref,icont_ref,co)
740 if(me.eq.king.or..not.out1file)
741 & write(iout,*)'Adding sidechains'
745 if (iti.ne.10 .and. itype(i).ne.ntyp1) then
748 do while (fail.and.nsi.le.maxsi)
749 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
752 if(fail) write(iout,*)'Adding sidechain failed for res ',
753 & i,' after ',nsi,' trials'
758 if (indpdb.eq.0) then
759 C Read sequence if not taken from the pdb file.
761 c print *,'nres=',nres
762 if (iscode.gt.0) then
763 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
765 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
767 C Convert sequence to numeric code
769 itype(i)=rescode(i,sequence(i),iscode)
771 C Assign initial virtual bond lengths
777 vbld(i+nres)=dsc(iabs(itype(i)))
778 vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
779 c write (iout,*) "i",i," itype",itype(i),
780 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
784 c print '(20i4)',(itype(i),i=1,nres)
787 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
789 if (itype(i).eq.ntyp1) then
793 else if (iabs(itype(i+1)).ne.20) then
795 else if (iabs(itype(i)).ne.20) then
802 if(me.eq.king.or..not.out1file)then
803 write (iout,*) "ITEL"
805 write (iout,*) i,itype(i),itel(i)
807 print *,'Call Read_Bridge.'
810 C 8/13/98 Set limits to generating the dihedral angles
815 read (inp,*) ndih_constr
816 if (ndih_constr.gt.0) then
818 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
819 if(me.eq.king.or..not.out1file)then
821 & 'There are',ndih_constr,' constraints on phi angles.'
823 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
827 phi0(i)=deg2rad*phi0(i)
828 drange(i)=deg2rad*drange(i)
830 if(me.eq.king.or..not.out1file)
831 & write (iout,*) 'FTORS',ftors
834 phibound(1,ii) = phi0(i)-drange(i)
835 phibound(2,ii) = phi0(i)+drange(i)
842 write (iout,'(a)') 'Boundaries in phi angle sampling:'
844 write (iout,'(a3,i5,2f10.1)')
845 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
851 cd print *,'NNT=',NNT,' NCT=',NCT
852 if (itype(1).eq.ntyp1) nnt=2
853 if (itype(nres).eq.ntyp1) nct=nct-1
855 if(me.eq.king.or..not.out1file)
856 & write (iout,'(a,i3)') 'nsup=',nsup
858 if (nsup.le.(nct-nnt+1)) then
859 do i=0,nct-nnt+1-nsup
860 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
866 & 'Error - sequences to be superposed do not match.'
869 do i=0,nsup-(nct-nnt+1)
870 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
872 nstart_sup=nstart_sup+i
878 & 'Error - sequences to be superposed do not match.'
881 if (nsup.eq.0) nsup=nct-nnt
882 if (nstart_sup.eq.0) nstart_sup=nnt
883 if (nstart_seq.eq.0) nstart_seq=nnt
884 if(me.eq.king.or..not.out1file)
885 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
886 & ' nstart_seq=',nstart_seq
888 c--- Zscore rms -------
889 if (nz_start.eq.0) nz_start=nnt
890 if (nz_end.eq.0 .and. nsup.gt.0) then
892 else if (nz_end.eq.0) then
895 if(me.eq.king.or..not.out1file)then
896 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
897 write (iout,*) 'IZ_SC=',iz_sc
899 c----------------------
902 if (.not.pdbref) then
903 call read_angles(inp,*38)
905 38 write (iout,'(a)') 'Error reading reference structure.'
907 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
908 stop 'Error reading reference structure'
912 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
922 call contact(.true.,ncont_ref,icont_ref,co)
924 c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
926 if (constr_dist.gt.0) call read_dist_constr
927 write (iout,*) "After read_dist_constr nhpb",nhpb
929 if(me.eq.king.or..not.out1file)
930 & write (iout,*) 'Contact order:',co
932 if(me.eq.king.or..not.out1file)
933 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
936 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
938 if(me.eq.king.or..not.out1file)
939 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
940 & icont_ref(1,i),' ',
941 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
945 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
946 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
947 & modecalc.ne.10) then
948 C If input structure hasn't been supplied from the PDB file read or generate
950 if (iranconf.eq.0 .and. .not. extconf) then
951 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
952 & write (iout,'(a)') 'Initial geometry will be read in.'
954 read(inp,'(8f10.5)',end=36,err=36)
955 & ((c(l,k),l=1,3),k=1,nres),
956 & ((c(l,k+nres),l=1,3),k=nnt,nct)
957 write (iout,*) "Exit READ_CART"
958 write (iout,'(8f10.5)')
959 & ((c(l,k),l=1,3),k=1,nres),
960 & ((c(l,k+nres),l=1,3),k=nnt,nct)
961 call int_from_cart1(.true.)
962 write (iout,*) "Finish INT_TO_CART"
965 dc(j,i)=c(j,i+1)-c(j,i)
966 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
970 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
972 dc(j,i+nres)=c(j,i+nres)-c(j,i)
973 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
979 call read_angles(inp,*36)
982 36 write (iout,'(a)') 'Error reading angle file.'
984 call mpi_finalize( MPI_COMM_WORLD,IERR )
986 stop 'Error reading angle file.'
988 else if (extconf) then
989 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
990 & write (iout,'(a)') 'Extended chain initial geometry.'
992 theta(i)=90d0*deg2rad
998 alph(i)=110d0*deg2rad
1001 omeg(i)=-120d0*deg2rad
1002 if (itype(i).le.0) omeg(i)=-omeg(i)
1005 if(me.eq.king.or..not.out1file)
1006 & write (iout,'(a)') 'Random-generated initial geometry.'
1010 if (me.eq.king .or. fg_rank.eq.0 .and. (
1011 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1015 call gen_rand_conf(itmp,*30)
1017 30 write (iout,*) 'Failed to generate random conformation',
1018 & ', itrial=',itrial
1019 write (*,*) 'Processor:',me,
1020 & ' Failed to generate random conformation',
1030 write (iout,'(a,i3,a)') 'Processor:',me,
1031 & ' error in generating random conformation.'
1032 write (*,'(a,i3,a)') 'Processor:',me,
1033 & ' error in generating random conformation.'
1036 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1041 & ' error in generating random conformation.'
1046 elseif (modecalc.eq.4) then
1047 read (inp,'(a)') intinname
1048 open (intin,file=intinname,status='old',err=333)
1049 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1050 & write (iout,'(a)') 'intinname',intinname
1051 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1053 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1055 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1057 stop 'Error opening angle file.'
1061 C Generate distance constraints, if the PDB structure is to be regularized.
1062 if (nthread.gt.0) then
1063 call read_threadbase
1066 if (me.eq.king .or. .not. out1file)
1068 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1069 write (iout,'(/a,i3,a)')
1070 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1071 write (iout,'(20i4)') (iss(i),i=1,ns)
1073 write(iout,*)"Running with dynamic disulfide-bond formation"
1075 write (iout,'(/a/)') 'Pre-formed links are:'
1081 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1082 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1088 if (ns.gt.0.and.dyn_ss) then
1092 forcon(i-nss)=forcon(i)
1099 dyn_ss_mask(iss(i))=.true.
1102 if (i2ndstr.gt.0) call secstrp2dihc
1103 c call geom_to_var(nvar,x)
1104 c call etotal(energia(0))
1105 c call enerprint(energia(0))
1106 c call briefout(0,etot)
1108 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1109 cd write (iout,'(a)') 'Variable list:'
1110 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1112 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1113 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1114 & 'Processor',myrank,': end reading molecular data.'
1118 c--------------------------------------------------------------------------
1119 logical function seq_comp(itypea,itypeb,length)
1121 integer length,itypea(length),itypeb(length)
1124 if (itypea(i).ne.itypeb(i)) then
1132 c-----------------------------------------------------------------------------
1133 subroutine read_bridge
1134 C Read information about disulfide bridges.
1135 implicit real*8 (a-h,o-z)
1136 include 'DIMENSIONS'
1140 include 'COMMON.IOUNITS'
1141 include 'COMMON.GEO'
1142 include 'COMMON.VAR'
1143 include 'COMMON.INTERACT'
1144 include 'COMMON.LOCAL'
1145 include 'COMMON.NAMES'
1146 include 'COMMON.CHAIN'
1147 include 'COMMON.FFIELD'
1148 include 'COMMON.SBRIDGE'
1149 include 'COMMON.HEADER'
1150 include 'COMMON.CONTROL'
1151 include 'COMMON.DBASE'
1152 include 'COMMON.THREAD'
1153 include 'COMMON.TIME1'
1154 include 'COMMON.SETUP'
1155 C Read bridging residues.
1156 read (inp,*) ns,(iss(i),i=1,ns)
1158 if(me.eq.king.or..not.out1file)
1159 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1160 C Check whether the specified bridging residues are cystines.
1162 if (itype(iss(i)).ne.1) then
1163 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1164 & 'Do you REALLY think that the residue ',
1165 & restyp(itype(iss(i))),i,
1166 & ' can form a disulfide bridge?!!!'
1167 write (*,'(2a,i3,a)')
1168 & 'Do you REALLY think that the residue ',
1169 & restyp(itype(iss(i))),i,
1170 & ' can form a disulfide bridge?!!!'
1172 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1177 C Read preformed bridges.
1179 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1181 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1184 C Check if the residues involved in bridges are in the specified list of
1185 C bridging residues.
1188 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1189 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1190 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1191 & ' contains residues present in other pairs.'
1192 write (*,'(a,i3,a)') 'Disulfide pair',i,
1193 & ' contains residues present in other pairs.'
1195 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1201 if (ihpb(i).eq.iss(j)) goto 10
1203 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1206 if (jhpb(i).eq.iss(j)) goto 20
1208 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1214 ihpb(i)=ihpb(i)+nres
1215 jhpb(i)=jhpb(i)+nres
1221 c----------------------------------------------------------------------------
1222 subroutine read_x(kanal,*)
1223 implicit real*8 (a-h,o-z)
1224 include 'DIMENSIONS'
1225 include 'COMMON.GEO'
1226 include 'COMMON.VAR'
1227 include 'COMMON.CHAIN'
1228 include 'COMMON.IOUNITS'
1229 include 'COMMON.CONTROL'
1230 include 'COMMON.LOCAL'
1231 include 'COMMON.INTERACT'
1232 c Read coordinates from input
1234 read(kanal,'(8f10.5)',end=10,err=10)
1235 & ((c(l,k),l=1,3),k=1,nres),
1236 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1239 c(j,2*nres)=c(j,nres)
1241 call int_from_cart1(.false.)
1244 dc(j,i)=c(j,i+1)-c(j,i)
1245 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1249 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1251 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1252 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1260 c----------------------------------------------------------------------------
1261 subroutine read_threadbase
1262 implicit real*8 (a-h,o-z)
1263 include 'DIMENSIONS'
1264 include 'COMMON.IOUNITS'
1265 include 'COMMON.GEO'
1266 include 'COMMON.VAR'
1267 include 'COMMON.INTERACT'
1268 include 'COMMON.LOCAL'
1269 include 'COMMON.NAMES'
1270 include 'COMMON.CHAIN'
1271 include 'COMMON.FFIELD'
1272 include 'COMMON.SBRIDGE'
1273 include 'COMMON.HEADER'
1274 include 'COMMON.CONTROL'
1275 include 'COMMON.DBASE'
1276 include 'COMMON.THREAD'
1277 include 'COMMON.TIME1'
1278 C Read pattern database for threading.
1279 read (icbase,*) nseq
1281 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1282 & nres_base(2,i),nres_base(3,i)
1283 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1285 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1286 c & nres_base(2,i),nres_base(3,i)
1287 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1291 if (weidis.eq.0.0D0) weidis=0.1D0
1300 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1301 write (iout,'(a,i5)') 'nexcl: ',nexcl
1302 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1305 c------------------------------------------------------------------------------
1306 subroutine setup_var
1307 implicit real*8 (a-h,o-z)
1308 include 'DIMENSIONS'
1309 include 'COMMON.IOUNITS'
1310 include 'COMMON.GEO'
1311 include 'COMMON.VAR'
1312 include 'COMMON.INTERACT'
1313 include 'COMMON.LOCAL'
1314 include 'COMMON.NAMES'
1315 include 'COMMON.CHAIN'
1316 include 'COMMON.FFIELD'
1317 include 'COMMON.SBRIDGE'
1318 include 'COMMON.HEADER'
1319 include 'COMMON.CONTROL'
1320 include 'COMMON.DBASE'
1321 include 'COMMON.THREAD'
1322 include 'COMMON.TIME1'
1323 C Set up variable list.
1329 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1331 ialph(i,1)=nvar+nside
1335 if (indphi.gt.0) then
1337 else if (indback.gt.0) then
1342 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1345 c----------------------------------------------------------------------------
1346 subroutine gen_dist_constr
1347 C Generate CA distance constraints.
1348 implicit real*8 (a-h,o-z)
1349 include 'DIMENSIONS'
1350 include 'COMMON.IOUNITS'
1351 include 'COMMON.GEO'
1352 include 'COMMON.VAR'
1353 include 'COMMON.INTERACT'
1354 include 'COMMON.LOCAL'
1355 include 'COMMON.NAMES'
1356 include 'COMMON.CHAIN'
1357 include 'COMMON.FFIELD'
1358 include 'COMMON.SBRIDGE'
1359 include 'COMMON.HEADER'
1360 include 'COMMON.CONTROL'
1361 include 'COMMON.DBASE'
1362 include 'COMMON.THREAD'
1363 include 'COMMON.TIME1'
1364 dimension itype_pdb(maxres)
1365 common /pizda/ itype_pdb
1367 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1368 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1369 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1371 do i=nstart_sup,nstart_sup+nsup-1
1372 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1373 cd & ' seq_pdb', restyp(itype_pdb(i))
1374 do j=i+2,nstart_sup+nsup-1
1376 ihpb(nhpb)=i+nstart_seq-nstart_sup
1377 jhpb(nhpb)=j+nstart_seq-nstart_sup
1379 dhpb(nhpb)=dist(i,j)
1382 cd write (iout,'(a)') 'Distance constraints:'
1387 cd if (ii.gt.nres) then
1392 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1393 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1394 cd & dhpb(i),forcon(i)
1398 c----------------------------------------------------------------------------
1400 implicit real*8 (a-h,o-z)
1401 include 'DIMENSIONS'
1402 include 'COMMON.MAP'
1403 include 'COMMON.IOUNITS'
1404 character*3 angid(4) /'THE','PHI','ALP','OME'/
1405 character*80 mapcard,ucase
1407 read (inp,'(a)') mapcard
1408 mapcard=ucase(mapcard)
1409 if (index(mapcard,'PHI').gt.0) then
1411 else if (index(mapcard,'THE').gt.0) then
1413 else if (index(mapcard,'ALP').gt.0) then
1415 else if (index(mapcard,'OME').gt.0) then
1418 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1419 stop 'Error - illegal variable spec in MAP card.'
1421 call readi (mapcard,'RES1',res1(imap),0)
1422 call readi (mapcard,'RES2',res2(imap),0)
1423 if (res1(imap).eq.0) then
1424 res1(imap)=res2(imap)
1425 else if (res2(imap).eq.0) then
1426 res2(imap)=res1(imap)
1428 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1430 & 'Error - illegal definition of variable group in MAP.'
1431 stop 'Error - illegal definition of variable group in MAP.'
1433 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1434 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1435 call readi(mapcard,'NSTEP',nstep(imap),0)
1436 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1438 & 'Illegal boundary and/or step size specification in MAP.'
1439 stop 'Illegal boundary and/or step size specification in MAP.'
1444 c----------------------------------------------------------------------------
1446 implicit real*8 (a-h,o-z)
1447 include 'DIMENSIONS'
1448 include 'COMMON.IOUNITS'
1449 include 'COMMON.GEO'
1450 include 'COMMON.CSA'
1451 include 'COMMON.BANK'
1452 include 'COMMON.CONTROL'
1454 character*620 mcmcard
1455 call card_concat(mcmcard)
1457 call readi(mcmcard,'NCONF',nconf,50)
1458 call readi(mcmcard,'NADD',nadd,0)
1459 call readi(mcmcard,'JSTART',jstart,1)
1460 call readi(mcmcard,'JEND',jend,1)
1461 call readi(mcmcard,'NSTMAX',nstmax,500000)
1462 call readi(mcmcard,'N0',n0,1)
1463 call readi(mcmcard,'N1',n1,6)
1464 call readi(mcmcard,'N2',n2,4)
1465 call readi(mcmcard,'N3',n3,0)
1466 call readi(mcmcard,'N4',n4,0)
1467 call readi(mcmcard,'N5',n5,0)
1468 call readi(mcmcard,'N6',n6,10)
1469 call readi(mcmcard,'N7',n7,0)
1470 call readi(mcmcard,'N8',n8,0)
1471 call readi(mcmcard,'N9',n9,0)
1472 call readi(mcmcard,'N14',n14,0)
1473 call readi(mcmcard,'N15',n15,0)
1474 call readi(mcmcard,'N16',n16,0)
1475 call readi(mcmcard,'N17',n17,0)
1476 call readi(mcmcard,'N18',n18,0)
1478 vdisulf=(index(mcmcard,'DYNSS').gt.0)
1480 call readi(mcmcard,'NDIFF',ndiff,2)
1481 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1482 call readi(mcmcard,'IS1',is1,1)
1483 call readi(mcmcard,'IS2',is2,8)
1484 call readi(mcmcard,'NRAN0',nran0,4)
1485 call readi(mcmcard,'NRAN1',nran1,2)
1486 call readi(mcmcard,'IRR',irr,1)
1487 call readi(mcmcard,'NSEED',nseed,20)
1488 call readi(mcmcard,'NTOTAL',ntotal,10000)
1489 call reada(mcmcard,'CUT1',cut1,2.0d0)
1490 call reada(mcmcard,'CUT2',cut2,5.0d0)
1491 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1492 call readi(mcmcard,'ICMAX',icmax,3)
1493 call readi(mcmcard,'IRESTART',irestart,0)
1494 c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1497 call reada(mcmcard,'DELE',dele,20.0d0)
1498 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1499 call readi(mcmcard,'IREF',iref,0)
1500 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1501 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1502 call readi(mcmcard,'NCONF_IN',nconf_in,0)
1503 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1504 write (iout,*) "NCONF_IN",nconf_in
1507 c----------------------------------------------------------------------------
1508 cfmc subroutine mcmfread
1509 cfmc implicit real*8 (a-h,o-z)
1510 cfmc include 'DIMENSIONS'
1511 cfmc include 'COMMON.MCMF'
1512 cfmc include 'COMMON.IOUNITS'
1513 cfmc include 'COMMON.GEO'
1514 cfmc character*80 ucase
1515 cfmc character*620 mcmcard
1516 cfmc call card_concat(mcmcard)
1518 cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
1519 cfmc write(iout,*)'MAXRANT=',maxrant
1520 cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1521 cfmc write(iout,*)'MAXFAM=',maxfam
1522 cfmc call readi(mcmcard,'NNET1',nnet1,5)
1523 cfmc write(iout,*)'NNET1=',nnet1
1524 cfmc call readi(mcmcard,'NNET2',nnet2,4)
1525 cfmc write(iout,*)'NNET2=',nnet2
1526 cfmc call readi(mcmcard,'NNET3',nnet3,4)
1527 cfmc write(iout,*)'NNET3=',nnet3
1528 cfmc call readi(mcmcard,'ILASTT',ilastt,0)
1529 cfmc write(iout,*)'ILASTT=',ilastt
1530 cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1531 cfmc write(iout,*)'MAXSTR=',maxstr
1532 cfmc maxstr_f=maxstr/maxfam
1533 cfmc write(iout,*)'MAXSTR_F=',maxstr_f
1534 cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
1535 cfmc write(iout,*)'NMCMF=',nmcmf
1536 cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1537 cfmc write(iout,*)'IFOCUS=',ifocus
1538 cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1539 cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
1540 cfmc call readi(mcmcard,'INTPRT',intprt,1000)
1541 cfmc write(iout,*)'INTPRT=',intprt
1542 cfmc call readi(mcmcard,'IPRT',iprt,100)
1543 cfmc write(iout,*)'IPRT=',iprt
1544 cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
1545 cfmc write(iout,*)'IMAXTR=',imaxtr
1546 cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
1547 cfmc write(iout,*)'MAXEVEN=',maxeven
1548 cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1549 cfmc write(iout,*)'MAXEVEN1=',maxeven1
1550 cfmc call readi(mcmcard,'INIMIN',inimin,200)
1551 cfmc write(iout,*)'INIMIN=',inimin
1552 cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1553 cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
1554 cfmc call readi(mcmcard,'NTHREAD',nthread,5)
1555 cfmc write(iout,*)'NTHREAD=',nthread
1556 cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1557 cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1558 cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
1559 cfmc write(iout,*)'MAXPERT=',maxpert
1560 cfmc call readi(mcmcard,'IRMSD',irmsd,1)
1561 cfmc write(iout,*)'IRMSD=',irmsd
1562 cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1563 cfmc write(iout,*)'DENEMIN=',denemin
1564 cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1565 cfmc write(iout,*)'RCUT1S=',rcut1s
1566 cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1567 cfmc write(iout,*)'RCUT1E=',rcut1e
1568 cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1569 cfmc write(iout,*)'RCUT2S=',rcut2s
1570 cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1571 cfmc write(iout,*)'RCUT2E=',rcut2e
1572 cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1573 cfmc write(iout,*)'DPERT1=',d_pert1
1574 cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1575 cfmc write(iout,*)'DPERT1A=',d_pert1a
1576 cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1577 cfmc write(iout,*)'DPERT2=',d_pert2
1578 cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1579 cfmc write(iout,*)'DPERT2A=',d_pert2a
1580 cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1581 cfmc write(iout,*)'DPERT2B=',d_pert2b
1582 cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1583 cfmc write(iout,*)'DPERT2C=',d_pert2c
1584 cfmc d_pert1=deg2rad*d_pert1
1585 cfmc d_pert1a=deg2rad*d_pert1a
1586 cfmc d_pert2=deg2rad*d_pert2
1587 cfmc d_pert2a=deg2rad*d_pert2a
1588 cfmc d_pert2b=deg2rad*d_pert2b
1589 cfmc d_pert2c=deg2rad*d_pert2c
1590 cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1591 cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
1592 cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1593 cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
1594 cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1595 cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1596 cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1597 cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1598 cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1599 cfmc write(iout,*)'RCUTINI=',rcutini
1600 cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
1601 cfmc write(iout,*)'GRAT=',grat
1602 cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1603 cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
1607 c----------------------------------------------------------------------------
1609 implicit real*8 (a-h,o-z)
1610 include 'DIMENSIONS'
1611 include 'COMMON.MCM'
1612 include 'COMMON.MCE'
1613 include 'COMMON.IOUNITS'
1615 character*320 mcmcard
1616 call card_concat(mcmcard)
1617 call readi(mcmcard,'MAXACC',maxacc,100)
1618 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1619 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1620 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1621 call readi(mcmcard,'MAXREPM',maxrepm,200)
1622 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1623 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1624 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1625 call reada(mcmcard,'E_UP',e_up,5.0D0)
1626 call reada(mcmcard,'DELTE',delte,0.1D0)
1627 call readi(mcmcard,'NSWEEP',nsweep,5)
1628 call readi(mcmcard,'NSTEPH',nsteph,0)
1629 call readi(mcmcard,'NSTEPC',nstepc,0)
1630 call reada(mcmcard,'TMIN',tmin,298.0D0)
1631 call reada(mcmcard,'TMAX',tmax,298.0D0)
1632 call readi(mcmcard,'NWINDOW',nwindow,0)
1633 call readi(mcmcard,'PRINT_MC',print_mc,0)
1634 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1635 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1636 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1637 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1638 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1639 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1640 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1641 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1642 if (nwindow.gt.0) then
1643 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1645 winlen(i)=winend(i)-winstart(i)+1
1648 if (tmax.lt.tmin) tmax=tmin
1649 if (tmax.eq.tmin) then
1653 if (nstepc.gt.0 .and. nsteph.gt.0) then
1654 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1655 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1657 C Probabilities of different move types
1658 sumpro_type(0)=0.0D0
1659 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1660 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1661 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1662 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1663 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1664 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1665 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1667 print *,'i',i,' sumprotype',sumpro_type(i)
1668 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1669 print *,'i',i,' sumprotype',sumpro_type(i)
1673 c----------------------------------------------------------------------------
1674 subroutine read_minim
1675 implicit real*8 (a-h,o-z)
1676 include 'DIMENSIONS'
1677 include 'COMMON.MINIM'
1678 include 'COMMON.IOUNITS'
1680 character*320 minimcard
1681 call card_concat(minimcard)
1682 call readi(minimcard,'MAXMIN',maxmin,2000)
1683 call readi(minimcard,'MAXFUN',maxfun,5000)
1684 call readi(minimcard,'MINMIN',minmin,maxmin)
1685 call readi(minimcard,'MINFUN',minfun,maxmin)
1686 call reada(minimcard,'TOLF',tolf,1.0D-2)
1687 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1688 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1689 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1690 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1691 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1692 & 'Options in energy minimization:'
1693 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1694 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1695 & 'MinMin:',MinMin,' MinFun:',MinFun,
1696 & ' TolF:',TolF,' RTolF:',RTolF
1699 c----------------------------------------------------------------------------
1700 subroutine read_angles(kanal,*)
1701 implicit real*8 (a-h,o-z)
1702 include 'DIMENSIONS'
1703 include 'COMMON.GEO'
1704 include 'COMMON.VAR'
1705 include 'COMMON.CHAIN'
1706 include 'COMMON.IOUNITS'
1707 include 'COMMON.CONTROL'
1708 c Read angles from input
1710 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1711 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1712 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1713 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1716 c 9/7/01 avoid 180 deg valence angle
1717 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1719 theta(i)=deg2rad*theta(i)
1720 phi(i)=deg2rad*phi(i)
1721 alph(i)=deg2rad*alph(i)
1722 omeg(i)=deg2rad*omeg(i)
1727 c----------------------------------------------------------------------------
1728 subroutine reada(rekord,lancuch,wartosc,default)
1730 character*(*) rekord,lancuch
1731 double precision wartosc,default
1734 iread=index(rekord,lancuch)
1735 if (iread.eq.0) then
1739 iread=iread+ilen(lancuch)+1
1740 read (rekord(iread:),*,err=10,end=10) wartosc
1745 c----------------------------------------------------------------------------
1746 subroutine readi(rekord,lancuch,wartosc,default)
1748 character*(*) rekord,lancuch
1749 integer wartosc,default
1752 iread=index(rekord,lancuch)
1753 if (iread.eq.0) then
1757 iread=iread+ilen(lancuch)+1
1758 read (rekord(iread:),*,err=10,end=10) wartosc
1763 c----------------------------------------------------------------------------
1764 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1767 integer tablica(dim),default
1768 character*(*) rekord,lancuch
1775 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1776 if (iread.eq.0) return
1777 iread=iread+ilen(lancuch)+1
1778 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1781 c----------------------------------------------------------------------------
1782 subroutine multreada(rekord,lancuch,tablica,dim,default)
1785 double precision tablica(dim),default
1786 character*(*) rekord,lancuch
1793 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1794 if (iread.eq.0) return
1795 iread=iread+ilen(lancuch)+1
1796 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1799 c----------------------------------------------------------------------------
1800 subroutine openunits
1801 implicit real*8 (a-h,o-z)
1802 include 'DIMENSIONS'
1805 character*16 form,nodename
1808 include 'COMMON.SETUP'
1809 include 'COMMON.IOUNITS'
1811 include 'COMMON.CONTROL'
1812 integer lenpre,lenpot,ilen,lentmp
1814 character*3 out1file_text,ucase
1817 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1818 call getenv_loc("PREFIX",prefix)
1820 call getenv_loc("POT",pot)
1821 call getenv_loc("DIRTMP",tmpdir)
1822 call getenv_loc("CURDIR",curdir)
1823 call getenv_loc("OUT1FILE",out1file_text)
1824 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1825 out1file_text=ucase(out1file_text)
1826 if (out1file_text(1:1).eq."Y") then
1829 out1file=fg_rank.gt.0
1834 if (lentmp.gt.0) then
1835 write (*,'(80(1h!))')
1836 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
1837 write (*,'(80(1h!))')
1838 write (*,*)"All output files will be on node /tmp directory."
1840 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1841 if (me.eq.king) then
1842 write (*,*) "The master node is ",nodename
1843 else if (fg_rank.eq.0) then
1844 write (*,*) "I am the CG slave node ",nodename
1846 write (*,*) "I am the FG slave node ",nodename
1849 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1850 lenpre = lentmp+lenpre+1
1852 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
1853 C Get the names and open the input files
1854 #if defined(WINIFL) || defined(WINPGI)
1855 open(1,file=pref_orig(:ilen(pref_orig))//
1856 & '.inp',status='old',readonly,shared)
1857 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1858 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1859 C Get parameter filenames and open the parameter files.
1860 call getenv_loc('BONDPAR',bondname)
1861 open (ibond,file=bondname,status='old',readonly,shared)
1862 call getenv_loc('THETPAR',thetname)
1863 open (ithep,file=thetname,status='old',readonly,shared)
1864 call getenv_loc('ROTPAR',rotname)
1865 open (irotam,file=rotname,status='old',readonly,shared)
1866 call getenv_loc('TORPAR',torname)
1867 open (itorp,file=torname,status='old',readonly,shared)
1868 call getenv_loc('TORDPAR',tordname)
1869 open (itordp,file=tordname,status='old',readonly,shared)
1870 call getenv_loc('FOURIER',fouriername)
1871 open (ifourier,file=fouriername,status='old',readonly,shared)
1872 call getenv_loc('ELEPAR',elename)
1873 open (ielep,file=elename,status='old',readonly,shared)
1874 call getenv_loc('SIDEPAR',sidename)
1875 open (isidep,file=sidename,status='old',readonly,shared)
1876 #elif (defined CRAY) || (defined AIX)
1877 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1879 c print *,"Processor",myrank," opened file 1"
1880 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1881 c print *,"Processor",myrank," opened file 9"
1882 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1883 C Get parameter filenames and open the parameter files.
1884 call getenv_loc('BONDPAR',bondname)
1885 open (ibond,file=bondname,status='old',action='read')
1886 c print *,"Processor",myrank," opened file IBOND"
1887 call getenv_loc('THETPAR',thetname)
1888 open (ithep,file=thetname,status='old',action='read')
1889 c print *,"Processor",myrank," opened file ITHEP"
1890 call getenv_loc('ROTPAR',rotname)
1891 open (irotam,file=rotname,status='old',action='read')
1892 c print *,"Processor",myrank," opened file IROTAM"
1893 call getenv_loc('TORPAR',torname)
1894 open (itorp,file=torname,status='old',action='read')
1895 c print *,"Processor",myrank," opened file ITORP"
1896 call getenv_loc('TORDPAR',tordname)
1897 open (itordp,file=tordname,status='old',action='read')
1898 c print *,"Processor",myrank," opened file ITORDP"
1899 call getenv_loc('SCCORPAR',sccorname)
1900 open (isccor,file=sccorname,status='old',action='read')
1901 c print *,"Processor",myrank," opened file ISCCOR"
1902 call getenv_loc('FOURIER',fouriername)
1903 open (ifourier,file=fouriername,status='old',action='read')
1904 c print *,"Processor",myrank," opened file IFOURIER"
1905 call getenv_loc('ELEPAR',elename)
1906 open (ielep,file=elename,status='old',action='read')
1907 c print *,"Processor",myrank," opened file IELEP"
1908 call getenv_loc('SIDEPAR',sidename)
1909 open (isidep,file=sidename,status='old',action='read')
1910 c print *,"Processor",myrank," opened file ISIDEP"
1911 c print *,"Processor",myrank," opened parameter files"
1913 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
1914 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1915 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1916 C Get parameter filenames and open the parameter files.
1917 call getenv_loc('BONDPAR',bondname)
1918 open (ibond,file=bondname,status='old')
1919 call getenv_loc('THETPAR',thetname)
1920 open (ithep,file=thetname,status='old')
1921 call getenv_loc('ROTPAR',rotname)
1922 open (irotam,file=rotname,status='old')
1923 call getenv_loc('TORPAR',torname)
1924 open (itorp,file=torname,status='old')
1925 call getenv_loc('TORDPAR',tordname)
1926 open (itordp,file=tordname,status='old')
1927 call getenv_loc('SCCORPAR',sccorname)
1928 open (isccor,file=sccorname,status='old')
1929 call getenv_loc('FOURIER',fouriername)
1930 open (ifourier,file=fouriername,status='old')
1931 call getenv_loc('ELEPAR',elename)
1932 open (ielep,file=elename,status='old')
1933 call getenv_loc('SIDEPAR',sidename)
1934 open (isidep,file=sidename,status='old')
1936 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1938 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1939 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1940 C Get parameter filenames and open the parameter files.
1941 call getenv_loc('BONDPAR',bondname)
1942 open (ibond,file=bondname,status='old',readonly)
1943 call getenv_loc('THETPAR',thetname)
1944 open (ithep,file=thetname,status='old',readonly)
1945 call getenv_loc('ROTPAR',rotname)
1946 open (irotam,file=rotname,status='old',readonly)
1947 call getenv_loc('TORPAR',torname)
1948 open (itorp,file=torname,status='old',readonly)
1949 call getenv_loc('TORDPAR',tordname)
1950 open (itordp,file=tordname,status='old',readonly)
1951 call getenv_loc('SCCORPAR',sccorname)
1952 open (isccor,file=sccorname,status='old',readonly)
1954 call getenv_loc('THETPARPDB',thetname_pdb)
1955 print *,"thetname_pdb ",thetname_pdb
1956 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
1957 print *,ithep_pdb," opened"
1959 call getenv_loc('FOURIER',fouriername)
1960 open (ifourier,file=fouriername,status='old',readonly)
1961 call getenv_loc('ELEPAR',elename)
1962 open (ielep,file=elename,status='old',readonly)
1963 call getenv_loc('SIDEPAR',sidename)
1964 open (isidep,file=sidename,status='old',readonly)
1966 call getenv_loc('ROTPARPDB',rotname_pdb)
1967 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
1972 C 8/9/01 In the newest version SCp interaction constants are read from a file
1973 C Use -DOLDSCP to use hard-coded constants instead.
1975 call getenv_loc('SCPPAR',scpname)
1976 #if defined(WINIFL) || defined(WINPGI)
1977 open (iscpp,file=scpname,status='old',readonly,shared)
1978 #elif (defined CRAY) || (defined AIX)
1979 open (iscpp,file=scpname,status='old',action='read')
1981 open (iscpp,file=scpname,status='old')
1983 open (iscpp,file=scpname,status='old',readonly)
1986 call getenv_loc('PATTERN',patname)
1987 #if defined(WINIFL) || defined(WINPGI)
1988 open (icbase,file=patname,status='old',readonly,shared)
1989 #elif (defined CRAY) || (defined AIX)
1990 open (icbase,file=patname,status='old',action='read')
1992 open (icbase,file=patname,status='old')
1994 open (icbase,file=patname,status='old',readonly)
1997 C Open output file only for CG processes
1998 c print *,"Processor",myrank," fg_rank",fg_rank
1999 if (fg_rank.eq.0) then
2001 if (nodes.eq.1) then
2004 npos = dlog10(dfloat(nodes-1))+1
2006 if (npos.lt.3) npos=3
2007 write (liczba,'(i1)') npos
2008 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2010 write (liczba,form) me
2011 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2012 & liczba(:ilen(liczba))
2013 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2015 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2017 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2018 & liczba(:ilen(liczba))//'.mol2'
2019 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2020 & liczba(:ilen(liczba))//'.stat'
2022 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2023 & //liczba(:ilen(liczba))//'.stat')
2024 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2027 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2028 & liczba(:ilen(liczba))//'.const'
2033 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2034 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2035 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2036 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2037 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2039 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2041 rest2name=prefix(:ilen(prefix))//'.rst'
2043 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2046 #if defined(AIX) || defined(PGI)
2047 if (me.eq.king .or. .not. out1file)
2048 & open(iout,file=outname,status='unknown')
2050 if (fg_rank.gt.0) then
2051 write (liczba,'(i3.3)') myrank/nfgtasks
2052 write (ll,'(bz,i3.3)') fg_rank
2053 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2058 open(igeom,file=intname,status='unknown',position='append')
2059 open(ipdb,file=pdbname,status='unknown')
2060 open(imol2,file=mol2name,status='unknown')
2061 open(istat,file=statname,status='unknown',position='append')
2063 c1out open(iout,file=outname,status='unknown')
2066 if (me.eq.king .or. .not.out1file)
2067 & open(iout,file=outname,status='unknown')
2069 if (fg_rank.gt.0) then
2070 write (liczba,'(i3.3)') myrank/nfgtasks
2071 write (ll,'(bz,i3.3)') fg_rank
2072 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2077 open(igeom,file=intname,status='unknown',access='append')
2078 open(ipdb,file=pdbname,status='unknown')
2079 open(imol2,file=mol2name,status='unknown')
2080 open(istat,file=statname,status='unknown',access='append')
2082 c1out open(iout,file=outname,status='unknown')
2085 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2086 csa_seed=prefix(:lenpre)//'.CSA.seed'
2087 csa_history=prefix(:lenpre)//'.CSA.history'
2088 csa_bank=prefix(:lenpre)//'.CSA.bank'
2089 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2090 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2091 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2092 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2093 csa_int=prefix(:lenpre)//'.int'
2094 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2095 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2096 csa_in=prefix(:lenpre)//'.CSA.in'
2097 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2100 write (iout,'(80(1h-))')
2101 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2102 write (iout,'(80(1h-))')
2103 write (iout,*) "Input file : ",
2104 & pref_orig(:ilen(pref_orig))//'.inp'
2105 write (iout,*) "Output file : ",
2106 & outname(:ilen(outname))
2108 write (iout,*) "Sidechain potential file : ",
2109 & sidename(:ilen(sidename))
2111 write (iout,*) "SCp potential file : ",
2112 & scpname(:ilen(scpname))
2114 write (iout,*) "Electrostatic potential file : ",
2115 & elename(:ilen(elename))
2116 write (iout,*) "Cumulant coefficient file : ",
2117 & fouriername(:ilen(fouriername))
2118 write (iout,*) "Torsional parameter file : ",
2119 & torname(:ilen(torname))
2120 write (iout,*) "Double torsional parameter file : ",
2121 & tordname(:ilen(tordname))
2122 write (iout,*) "SCCOR parameter file : ",
2123 & sccorname(:ilen(sccorname))
2124 write (iout,*) "Bond & inertia constant file : ",
2125 & bondname(:ilen(bondname))
2126 write (iout,*) "Bending parameter file : ",
2127 & thetname(:ilen(thetname))
2128 write (iout,*) "Rotamer parameter file : ",
2129 & rotname(:ilen(rotname))
2130 write (iout,*) "Thetpdb parameter file : ",
2131 & thetname_pdb(:ilen(thetname_pdb))
2132 write (iout,*) "Threading database : ",
2133 & patname(:ilen(patname))
2135 &write (iout,*)" DIRTMP : ",
2137 write (iout,'(80(1h-))')
2141 c----------------------------------------------------------------------------
2142 subroutine card_concat(card)
2143 implicit real*8 (a-h,o-z)
2144 include 'DIMENSIONS'
2145 include 'COMMON.IOUNITS'
2147 character*80 karta,ucase
2149 read (inp,'(a)') karta
2152 do while (karta(80:80).eq.'&')
2153 card=card(:ilen(card)+1)//karta(:79)
2154 read (inp,'(a)') karta
2157 card=card(:ilen(card)+1)//karta
2160 c----------------------------------------------------------------------------------
2162 implicit real*8 (a-h,o-z)
2163 include 'DIMENSIONS'
2164 include 'COMMON.CHAIN'
2165 include 'COMMON.IOUNITS'
2167 open(irest2,file=rest2name,status='unknown')
2168 read(irest2,*) totT,EK,potE,totE,t_bath
2170 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2173 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2176 read (irest2,*) iset
2181 c---------------------------------------------------------------------------------
2182 subroutine read_fragments
2183 implicit real*8 (a-h,o-z)
2184 include 'DIMENSIONS'
2188 include 'COMMON.SETUP'
2189 include 'COMMON.CHAIN'
2190 include 'COMMON.IOUNITS'
2192 include 'COMMON.CONTROL'
2193 read(inp,*) nset,nfrag,npair,nfrag_back
2194 if(me.eq.king.or..not.out1file)
2195 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2196 & " nfrag_back",nfrag_back
2198 read(inp,*) mset(iset)
2200 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2202 if(me.eq.king.or..not.out1file)
2203 & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2204 & ifrag(2,i,iset), qinfrag(i,iset)
2207 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2209 if(me.eq.king.or..not.out1file)
2210 & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2211 & ipair(2,i,iset), qinpair(i,iset)
2214 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2215 & wfrag_back(3,i,iset),
2216 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2217 if(me.eq.king.or..not.out1file)
2218 & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2219 & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2224 c-------------------------------------------------------------------------------
2225 subroutine read_dist_constr
2226 implicit real*8 (a-h,o-z)
2227 include 'DIMENSIONS'
2231 include 'COMMON.SETUP'
2232 include 'COMMON.CONTROL'
2233 include 'COMMON.CHAIN'
2234 include 'COMMON.IOUNITS'
2235 include 'COMMON.SBRIDGE'
2236 integer ifrag_(2,100),ipair_(2,100)
2237 double precision wfrag_(100),wpair_(100)
2238 character*500 controlcard
2239 c write (iout,*) "Calling read_dist_constr"
2240 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2242 call card_concat(controlcard)
2243 call readi(controlcard,"NFRAG",nfrag_,0)
2244 call readi(controlcard,"NPAIR",npair_,0)
2245 call readi(controlcard,"NDIST",ndist_,0)
2246 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2247 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2248 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2249 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2250 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2251 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2252 c write (iout,*) "IFRAG"
2254 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2256 c write (iout,*) "IPAIR"
2258 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2262 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2263 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2264 & ifrag_(2,i)=nstart_sup+nsup-1
2265 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2267 if (wfrag_(i).gt.0.0d0) then
2268 do j=ifrag_(1,i),ifrag_(2,i)-1
2269 do k=j+1,ifrag_(2,i)
2270 c write (iout,*) "j",j," k",k
2272 if (constr_dist.eq.1) then
2277 forcon(nhpb)=wfrag_(i)
2278 else if (constr_dist.eq.2) then
2279 if (ddjk.le.dist_cut) then
2284 forcon(nhpb)=wfrag_(i)
2291 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2294 if (.not.out1file .or. me.eq.king)
2295 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2296 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2298 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2299 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2306 if (wpair_(i).gt.0.0d0) then
2314 do j=ifrag_(1,ii),ifrag_(2,ii)
2315 do k=ifrag_(1,jj),ifrag_(2,jj)
2319 forcon(nhpb)=wpair_(i)
2320 dhpb(nhpb)=dist(j,k)
2322 if (.not.out1file .or. me.eq.king)
2323 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2324 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2326 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2327 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2334 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2335 if (forcon(nhpb+1).gt.0.0d0) then
2337 dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
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 c-------------------------------------------------------------------------------
2353 subroutine flush(iu)
2358 subroutine flush(iu)
2363 c------------------------------------------------------------------------------
2364 subroutine copy_to_tmp(source)
2365 include "DIMENSIONS"
2366 include "COMMON.IOUNITS"
2367 character*(*) source
2368 character* 256 tmpfile
2372 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
2373 inquire(file=tmpfile,exist=ex)
2375 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
2376 & " to temporary directory..."
2377 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
2378 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
2382 c------------------------------------------------------------------------------
2383 subroutine move_from_tmp(source)
2384 include "DIMENSIONS"
2385 include "COMMON.IOUNITS"
2386 character*(*) source
2389 write (*,*) "Moving ",source(:ilen(source)),
2390 & " from temporary directory to working directory"
2391 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
2392 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
2395 c------------------------------------------------------------------------------
2396 subroutine random_init(seed)
2398 C Initialize random number generator
2400 implicit real*8 (a-h,o-z)
2401 include 'DIMENSIONS'
2404 logical OKRandom, prng_restart
2406 integer iseed_array(4)
2408 include 'COMMON.IOUNITS'
2409 include 'COMMON.TIME1'
2410 include 'COMMON.THREAD'
2411 include 'COMMON.SBRIDGE'
2412 include 'COMMON.CONTROL'
2413 include 'COMMON.MCM'
2414 include 'COMMON.MAP'
2415 include 'COMMON.HEADER'
2416 include 'COMMON.CSA'
2417 include 'COMMON.CHAIN'
2418 include 'COMMON.MUCA'
2420 include 'COMMON.FFIELD'
2421 include 'COMMON.SETUP'
2422 iseed=-dint(dabs(seed))
2423 if (iseed.eq.0) then
2424 write (iout,'(/80(1h*)/20x,a/80(1h*))')
2425 & 'Random seed undefined. The program will stop.'
2426 write (*,'(/80(1h*)/20x,a/80(1h*))')
2427 & 'Random seed undefined. The program will stop.'
2429 call mpi_finalize(mpi_comm_world,ierr)
2431 stop 'Bad random seed.'
2434 if (fg_rank.eq.0) then
2438 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
2439 OKRandom = prng_restart(me,iseed)
2442 tmp=65536.0d0**(4-i)
2443 iseed_array(i) = dint(seed/tmp)
2444 seed=seed-iseed_array(i)*tmp
2447 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
2448 & (iseed_array(i),i=1,4)
2449 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
2450 & (iseed_array(i),i=1,4)
2451 OKRandom = prng_restart(me,iseed_array)
2454 c r1 = prng_next(me)
2455 r1=ran_number(0.0D0,1.0D0)
2457 & write (iout,*) 'ran_num',r1
2458 if (r1.lt.0.0d0) OKRandom=.false.
2460 if (.not.OKRandom) then
2461 write (iout,*) 'PRNG IS NOT WORKING!!!'
2462 print *,'PRNG IS NOT WORKING!!!'
2465 call mpi_abort(mpi_comm_world,error_msg,ierr)
2468 write (iout,*) 'too many processors for parallel prng'
2469 write (*,*) 'too many processors for parallel prng'
2477 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)